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ABSTRACT 



A two dimensional, incompressible viscous inviscid interaction computer code, de- 
signed to compute cascade flows, was investigated. Comparison of the flow character- 
istics predicted by the code with experimentally available data indicates that the code 
predicts reasonably well flow parameters on lightly loaded cascades. However, the code 
fails to predict correctly the actual boundary layer development and the velocity dis- 
tribution for highly loaded cascades. It is concluded that further improvement of the 
code is needed and recommendations are presented to achieve the required improve- 
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I. INTRODUCTION 



The need lor better and more efficient gas turbines requires the availability of cheap 
and reliable design tools for blades used in compressors and turbines. Computational 
methods are the preferred choice for such a design tool, considering the cost and com- 
plexity of wind tunnel experiments. 

Among the computational methods available today, the logical choice seems to be 
a computer code that can solve directly the full Xavier Stokes equations. However, given 
the state of the art in both algorithms and computer hardware, such Xavier Stokes 
solvers are restricted only to supercomputers, and even then the computation time is 
quite long. 

In order to enable fast and efficient computations, the viscous inviscid interaction 
code was developed by Cebeci [Ref. 1]. The approach used in this code is to solve the 
outer flow field using potential methods, and solving the boundary layer flow using a 
boundary layer method subject to an interaction law, that couples the inner and the 
outer flows. This interaction law is needed because classical boundary layer methods fail 
in areas of flow reversal and separation, which are very common in real life flows. 

The viscous inviscid interaction code was originally developed, and successfully used, 
for flows about single airfoils. It was later adapted to cascade flows. 

In this thesis the applicability of the code to cascades was investigated by comparing 
its results to experimental data. It was found that although the code can reasonably 
predict experimental results in some cases, it still needs improvements before it can be 
applied generally as a reliable design tool. 

A major restriction in improving the code is the lack of a wide data base of appro- 
priate experimental results. Some of the key elements in the code, like transition and 
turbulence modelling, are based on empirical correlations, and more detailed and accu- 
rate experiments should be performed, before a better understanding of these phenom- 
ena can be achieved. 

In the following, the theoretical background of the code is presented in Chapter II. 
a description of the code in Chapter III, the results and discussions are presented in 
Chapter IV and the conclusions and recommendations in Chapter V. A listing of the 
computer program is given in Appendix A. 



II. CASCADE FLOW PROBLEM FORMULATION 



This chapter outlines the theoretical background of the viscous inviscid interactive 
method used in the computer code to investigate cascade flows. Only the major steps in 
the mathematical developments will be described here. A detailed description of the 
theory and the numerical methods is given by Cebeci and Bradshaw [Ref. 2] and by 
Krainer [Ref. 3] on which this chapter is based. 

A. INVISCID FLOW METHOD 

Inviscid flow is the first building block of the flow and is solved using the panel 
method. The incompressible two dimensional outer flow must satisfy the Laplace 
equation: 

V 2 ci> = 0 , 

subject to the boundary conditions on the surface of the blade: 



where the commonlv used boundarv condition of zero normal velocitv on the surface is 
replaced by a specified blowing velocity v„ to allow for the effect of the boundary layer 
on the outer flow. 

In addition, the Kutta condition must be satisfied, in order to prevent the existence 
of discontinuous pressure distribution near the trailing edge of the blade. 

Since the Laplace equation is linear, a solution to a complex flow field can be built 
by superposition of solutions of elementary flows. The elementary flows used in the 
panel method are the uniform parallel flow and flows about singularities (sources and 
vortices). 

The panel method is based on replacing each blade by a distribution of sources and 
vortices on its surface. The surface is divided into a finite number of straight segments, 
called panels. 

On each panel, a uniform source distribution and a uniform vorticity distribution is 
assumed. The source strength at each panel is set to satisfy the boundary condition at 
the midpoint of the panel (called the control point), and so, in general the source 
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strength will vary from panel to panel. The vorticity strength is assumed to be the same 
for all the panels and is set to satisfy the Kutta condition. 

The cascade is defined as an infinite row of similar blades, each one modelled by 
panels of source and vortex distributions. The How at each point is found by summing 
the contributions of all the singularities on all the blades , and the uniform parallel flow. 

A useful concept in dealing with such flows is the concept of influence coefficients. 
An influence coefficient is defined as the velocity at a point induced by a unit strength 
singularity placed at some other field point. The influence coefficients are a function of 
geometry and so can be computed for a given cascade and a given choice of panel ge- 
ometry. 

Using the influence coefficients, the normal and tangential velocities at each control 
point can be written as a function of the unknown source strength of each panel and the 
unknown vortex strength. Using the boundary conditions, by equating the normal ve- 
locity at each control point to the prescribed blowing velocity v„, and using the Kutta 
condition (which requires equal velocities on the upper and on the lower panels at the 
trailing edge), a system of linear equations is constructed. 

By solving the system of equations, the strength of the sources and vortices is found, 
and the velocities (and the pressures) can be computed everywhere in the flow field. 

The velocity distribution on the suiface of the blade, computed by the panel method, 
is used as the boundary condition for the boundary layer flow calculations. 

It should be noted that in the panel method as used in the present computer code, 
there is no modelling of the wake, and its effect on the flow field is ignored. 

B. VISCOUS FLOW METHOD 

Viscous flow is the second building block of the cascade flow and it is applied to the 
thin boundary layer near the blade surface. 

1. Boundary Layer Theory 

The boundary layer concept, first suggested by L. Prandtl, assumes that the flow 
field can be divided into an outer flow where the viscous effects are negligible compared 
to inertia effects, and a thin layer close to the surface where the viscous effects cannot 
be neglected. The complete presentation of the boundary layer theory, and the devel- 
opment of the boundary layer equations, is given by Schlichting [Ref. 4). 

Under the assumptions of two dimensional incompressible thin boundary layer 
flow, the Xavier Stokes equations and the continuity equation reduce to: 
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= 0, 



CU _TT 

cx 3y 



cu , cu du e i c f , 8u 

u — I- v — = u e — h v - 7 — b ~ 

cx cy ax cy \ ov 



with the boundary conditions: 

y = 0 u(a\0) = 0 , v(x,0) = 0 , 

y — y'e u(xy e ) = u e ( x ) . 

v, 

where b denotes 1 + — . 

Writing the velocity components in terms of a stream function T : 

dip 

u = -r~ , 

cy 

8 ip 

cx 

This eliminates the continuity equation (which the stream function satisfies by defi- 
nition). 



Introducing the Falkner Skan transformation: 



/ u e 

} i v va' y ’ 



J[xy) = - ■ J=r 1 p(xy) , 



the momentum equation and the boundary conditions transform to: 

(*/")' + JR Y L ff" + 



V = 0 f(x, 0) = 0 = 0 , 
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// = >1 e 



f’{x, }J e ) = 1 , 



where m is defined by: 



m = 




<K 

dx 



The third order differential equation can be reduced to a system of first order 
differential equations by introduction of two new variables U and V : 

U = /' . 

V = L' , 



{b yy + + ,„(! - U 2 ) = ,(t f- - V -f ) , 



with the boundarv conditions: 



1 / = 0 U(.v,0) = 0 ,J[x, 0) = 0 , 



n = n e 






The next step is to use a finite difference approach to solve the equations. The 
box method is applied using central differencing in both the x and t] directions, and 
satisfying the equations midway between nodes. 

Applying the box method results in a system of nonlinear equations in the un- 
known variables (which are f, U and V in each node along the rj direction at the current 
x station). 

In order to solve the nonlinear system the Newton iterative procedure is used, 
linearizing the equations first about the solution at the adjacent upstream station, and 
then about the preceding iteration. The linearization is performed by letting: 




|T i,K t K— 1 



vp = vp- 1 



+#\ 

+ «Lp, 

+ 0Vp, 
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where: 



• / denotes location in the x direction 

• j denotes location in the y ( ?/ ) direction 

• k indicates the iteration counter 

This linearization results in a system of linear equations for the unknown in- 
crements: Sf / K , 5V,’;* and SV 1 ;* . 

This system of equations is solved repeatedly until the changes in the unknowns 
are small enough. Since the system is block tridiagonal, Keller's block elimination 
method is used. 

The method described so far, is a direct boundary' layer method. It can be used 
as long as the flow does not separate. Whenever separation or flow reversal occurs, and 
a zero skin friction coefficient is encountered, the equations become singular and the 
calculations will break down. 

2. Interactive Boundary Layer Method 

The interactive boundary layer method is designed to overcome the difficulties 
encountered at regions of flow reversal and separations. In such areas the external ve- 
locity is substantially changed by the viscous effects and can no longer be considered as 
a known boundary condition for the boundary layer flow. 

The general approach to the solution is the same as for the direct method but. 
since the outer flow is unknown, the velocity at the edge of the boundary layer is written 
as: 



"(-TJe) = u el + n 



d t x*v dl - 

~J7 (l<e d ) 7 

d£ e x — £ 



where: 

1. is total velocity at the edge of the boundary 7 layer. 

2. u fI (.v) is the velocity as computed by the inviscid method. 



/% 




is the Hilbert integral. 
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The numerical solution of the boundary layer equations follows the same steps 
as for the direct method, but with some changes. 

The transformation of the stream function and the y coordinate uses a constant 
velocity u 0 as a scaling factor, and a scaled velocity w is introduced: 

/~ 

V vjc y * 

A*, '/) = ■ — — Mxy), 

N /u 0 v.v 

u e( x -y) 

w = . 

Uq 

Using this transformation, the boundary layer equations become a system of 
first order differential equations: 



/' = U , 

U' = V , 
W' = 0 , 



{bV)'+4rJV + x\V-£r- 

2 cx 



= x\ 



d\J 

dx 



- V 



3L 

cx 



» 



with the boundary conditions: 

n = o u(x,0) = o ,j[x,o) = o , 
V = Ve U(*, Ve) = W (*> n e ) . 



w(x, i] e ) 







jc — <5 



The finite difference box method is used to solve the equations, in the same way 
as it was used for the direct case, but with two additions: 
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1. In areas of flow reversal the term ucu'8: x is omitted to assure stable integration 
(the FLARE approximation). 

2. The edge velocity, Wj (where J denotes the edge station) which involves integration, 
is approximated by : 



where g , and c„ are obtained from the numerical approximation to the Hilbert in- 
tegral (which will be presented in the next section ). 

By using central diffrencing to approximate the differential equations, a system 
of nonlinear algebraic equations is obtained for the unknown variables (which are 
fj, Uj , V; and W'). To solve the system of equations, the system is linearized by the 
Newton iterative procedure, and the resulting linear system is solved (for the new un- 
known variables which are the increments Sfy * , dl'y , and <5\V' ,r ). 

The solution of the system is repeated until the change in the increments is 
negligible compared to the preceding iteration, and the whole process is performed again 
at the next downstream station. 

3. Interactive Model 

The interactive model is used to couple the boundary layer to the external flow. 
It is needed in areas where strong interaction occurs, and both the boundary layer and 
the outer flow must be solved simultaneously. The interaction model provides the outer 
boundary condition to the boundary layer calculations by adding a correction term to 
the external velocity computed by the inviscid flow method. 

The external velocity is assumed to consist of a potential flow term ( u el (x ) ) and 
a correction term due to viscous effects ( u ti (.v) ): 



The viscous effect is obtained by a surface distribution of sources on the blade (a con- 
cept first suggested by Lighthill [Ref. 5]). The normal velocities at the surface of the 
blade, induced by these sources, displace the streamlines from the surface in the same 
way that the actual boundary layer displaces them: 



w; = gl + c„{w; , h ~fy , 



w„(.y) = «„(.<-) + uJX) . 



dS (x) 
dx 




Where v(.v, S’) is the normal velocity at the displaced surface. 



8 



Assuming that the surface can be approximated to be a flat plate, the normal 
velocity will be half the local source strength a (x). Assuming also that the inviscid ve- 
locity does not change across the boundary layer, the local source strength will be: 



ff(.v) 

? 



= v(*,0) = v{x, <5 ) 



rs 



"0 



cv 

cv 



dv = 






The local horizontal velocity induced by the source distribution, is the cor- 
rection term to the inviscid velocity, and can be represented by the Hilbert integral: 



_L 

71 



C x b 



<*(£) 



.v — 



dc = 



71 



d . d! 
— M ) 



The integration is carried out on all the sources on the surface, since the horizontal 
velocity is influenced by all the sources. 

The Hilbert integral is then approximated by a finite series: 



1 

71 





K 

c ik^ li e^ ) * 

/:=! 



Where c, k is a matrix of interaction coefficients which are functions of the geometry only 
( / denotes the chordwise position where u et) is evaluated and k is the location of the 
source which effects u e6 ). 

Since the computation of u e6 involves values of <5* downstream of the current x 
location, which are not known yet, these terms are taken from the previous iteration 
using a relaxation formula. 

4. Turbulence Model 

The turbulence model used here is the aleebraic eddv viscositv formulation of 

W 0* * 

Cebeci and Smith [Ref. 6]. According to the model used in the present computer code, 
the eddy viscosity v, is defined by two different expressions, for the inner region and for 
the outer region: 
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Where: 



v 



t 




0.4 y 



[■ 




cu 

oy 



y'n 



for 0<v<v c , 




- u)dyy tr y 



for y c <y <5 . 



A = 



26v 




1/2 






1 



1 + 


5.5(y!d) 6 


> 


0.0168 




1 - p 


du’8x 


2.5 


cu'cy 





P 



6 

1 + 2R t (2 - R t ) 



for R T < 1 , 



P = 



1 + R t 
R x 



for R t > 1 , 



Rt 



1 W 



( - Wv’) 



max 



The distance from the wall to the point between the two regions,^, is chosen such that 
the viscosity will be continuous. 
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The intermittency factor, y r , is defined by: 



Vir = 1 - exp 



u. 



Gy V 



RJi’ 34 (* 



- x tr ) 



dl 

IL 



Where: 

• R„, is the Reynolds number based on external velocity and transition location. 

• G v is an empirical constant, originally assigned the value 1200. 



Cebeci and Bradshaw [Ref. 2, p.246] described a different expression for the 
variable A in the inner region viscosity formula: 



A = 



26v 

0 - 

\ c y /max 



Where: 



vu e du e 



( , \ 3 / 2 ^ A * 
V oy J 



This version of the turbulence model was not implemented in the original computer 
code. During the work on this thesis, the effect of the modified turbulence model was 
investicated. 

A different intermittency distribution was implemented successfully by Rodi and 
Schonung [Ref. 7 ] for transition over separation bubbles. They used for G, the ex- 
pression: 



C = 100 

y exp(0.99 Tu) ' 

Where Tu is the turbulence level in the free flow. This intermittency model was also in- 
vestigated during the work on this thesis. 
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5. Transition 

The prediction of transition from laminar to turbulent llow is very difficult and 
has to rely on empirical correlations. The relation used here to predict the onset of 
transition is a combination of Michel's method and the e 9 method, and is given by 
Cebeci and Bradshaw [Ref. 2 , p. 153]: 



Where: 

1. is the Reynolds number based on the momentum thickness at the onset of 
transition. 

2. R, is the Revnolds number based on x at the onset of transition. 

In the computer code, if a laminar separation is detected before transition oc- 
curs, the onset of transition is assumed at the point of laminar separation. 
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III. DESCRIPTION OF THE COMPUTER CODE 



The computer code used here to investigate cascade flows was written by Cebeci, 
and is based on the numerical formulation that was outlined in the previous chapters. 
In this chapter the general structure and the major subroutines of the code will be de- 
scribed. 

A. GENERAL STRUCTURE OF THE MAIN PROGRAM 

The main program reads in the cascade data (blade coordinates, spacing and stagger 
angle), the flow data (inlet angle and Reynolds number), and transition parameters. The 
transition onset on each surface of the blade can be computed by the program, or can 
be input by the user. The intermittency parameter G should be specified by the user. 

The program then calls subroutine POTNL to compute the outer inviscid flow field 
for the first cycle. The output of subroutine POTNL is the external velocity distribution 
on the surface of the blades. This velocity distribution is then transferred to subroutine 
CASBLP, which calculate the boundary laver flow. 

Subroutine CASBLP returns the displacement thickness distribution and the blow- 
ing velocity distribution on the blades to the main program. This data is then transferred 
back to subroutine POTNL to the next cycle of calculations. 

The program repeats the cycles of calculations by calling the two subroutines, until 
the specified number of cycles is reached, or until a convergence criterion is satisfied. 

B. DESCRIPTION OF THE SUBROUTINES 
1. Subroutine POTNL 

This subroutine solves the inviscid outer flow by using the panel method. The 
subroutine calculates the influence coefficients and calculates the velocities subject to the 
boundary 7 conditions. 

The velocities are evaluated on the displaced surface (the surface created by 
adding the displacement thickness to the original surface of the blade). The input to this 
subroutine includes the cascade geometry 7 , the blowing velocity and the displacement 
thickness (for the first cycle both the displacement thickness and the blowing velocity- 
are taken to be zero). 
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2. Subroutine CASBLP 

This subroutine, called by the MAIN program, receives the blade geometry and 
the velocity distribution as input. 

It transforms the x.y blade coordinates to the chordwise tangential coordinates 
and smooths the velocity data (during the work on this thesis it was found that 
smoothing the velocity data prevents the detection of the separation bubble near the 
leading edge, and therefore it was eliminated). The subroutine then calls subroutine 
COMPBL for further calculations. 

3. Subroutine COMPBL 

This subroutine finds the stagnation point and controls the generation of the 
boundary’ layer calculation grid for each surface (the grid starts at the stagnation point 
and includes 91 points in the chordwise direction for the upper surface and 71 points on 
the lower surface). 

The subroutine then calls subroutine BL2D which calculates the boundary layer 
parameters for each surface (BL2D is called twice, first for the upper surface and then 
for the lower surface). 

4. Subroutine BL2D 

This subroutine computes the displacement thickness and the blowing velocity 
and returns them back to the calling subroutine (COMBL) in arrays compatible with the 
potential flow calculations (one array that contains all the points of the blade, first the 
lower surface starting at the trailing edge and proceeding forward, and then the upper 
surface, starting at the leading edge and proceeding backwards). 

BL2D calls the following subroutines: 

1. Subroutine INPUT which calculates the following: 

a. NS, the switching point between direct and interactive boundary layer calcu- 
lations (this point is set at the first pressure peak when the blade is scanned from 
leading edge towards the trailing edge) 

b. NTR, transition location (only if the transition location is an input. Otherwise 
it is calculated by subroutine TRNS). 

c. GMTR. the distribution of the intermittency factor y tr . 

In addition this subroutine generates the boundary’ layer grid in the t/ direction and 

the initial velocity profile, by calling subroutine IN'TL. 

2. CALCIJ, calculates the c„ coefficients used in the Hilbert integral approximation. 

3. EDDY, calculates the eddy viscosity (called only after transition has been de- 
tected). 
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4. COEFTR. calculates the coefiicients of the boundary layer finite difference 
equations in transformed form ( for the direct method calculations). 

5. SOLYE3, solves the linearized boundary layer equations for the F,L and V vari- 
ables by computing the increments OF. <)U and OV 

The subroutine then checks the convergence of the Newton iterations and re- 
peats the calculations if needed. If the subroutine detects flow separation or if it reaches 
the switching point NS, subroutine MA1N2 is called for the interactive method calcu- 
lations. Otherwise, the subroutine proceeds to the next chordwise point of the grid (NX) 
and repeats the calculations. 

5. Subroutine MAIN2 

This subroutine calculates the boundary layer parameters by the interactive 
method. The subroutine performs the following steps: 

1. It first calls subroutines JOINT and COMGI to compute the interaction coeffi- 
cients. 

2. In regions of laminar flow it calls the following subroutines: 

a. COEF. which calculates the coefficients of the boundary layer finite differences 
equations. 

b. SOLV4, solves for the variables F, U, V and W bv computing the increments 
5 F, <>U , <5V and d\V . 

c. TITANS, to check if the condition for transition is satisfied (it also checks for 
laminar separation and initiates transition at the point of laminar separation if 
it is detected). 

The subroutine then checks for convergence of Newton iterations and repeats the 
calculations as needed. 

3. In regions of turbulent flows the subroutine calls the following subroutines: 

a. EDDY, to compute the eddy viscosity parameter B (B = 1 + v f /v ). 

b. COEF and SOLV4, the same as for laminar flow. 



6. Subroutine OUTPUT 

This subroutine computes the boundary layer parameters. It is called with a 
parameter "INDEX" which determines the type of calculations: 

1. For INDEX=1 the computations relates to transformed coordinates (direct 
boundary layer method) using the relations: 



c f ~ 



2 V(1.2) B(l,2) 

v%7 
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<5* = -= (>;(A7>) ~ F (AT)) , 

V 1'or 



v 



NP 




Where V(l,2) and B( 1,2) are the velocity gradient and the viscosity parameter at the 
surface, respectively, and j/(NP) and F(NP) are ?/ and F evaluated at the edge of the 
boundary layer. 

2. For INDEX = 2 the subroutine calculates the boundary layer parameters for semi- 
-transformed coordinates (interactive boundary layer method) using the relations: 



For NX > NTR (after transition has been detected), subroutine SMPSON is 
called (subroutine SMPSON’ calculates the coefficient of the outer region eddy viscosity). 
The subroutine then prints out the velocity profiles at the required stations. 

7. Subroutine TRANS 

This subroutine calculates the transition location based on the Michel criterion 
or based on laminar separation (whichever occurs first). If transition has been detected 
the intermittency distribution is calculated for all the remaining points of the surface. 

8. Subroutine FILLUP 

This subroutine increases the number of points in the boundary layer grid (in 
the t] direction) as needed. It also fills up the arrays of F, U, B, W and V between the 
edge of the boundary layer to the end of the arrays (with V = 0, \V,B and U, with the last 
values that they had in the edge of the boundary layer and F as the integral of U). 



2 V(1.2) B( 1 .2) 



' f JZJK [W(A'P)] 2 ’ 
V(1 .2) 



u e = U(A T) , 

D = (U?7 - F) v \7 , 
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9. Subroutine EDDY 



This subroutine calculates the eddy viscosity using the Cebeci--Smith two layer 
eddy viscosity formula. It receives the vectors U,V and ?/ at a point and computes the 
viscosity vector B. 

10. Subroutine INTL 

This subroutine generates the boundary layer grid in the ?/ direction. It sets the 
number of grid points and generates the initial velocity profile. 
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IV. RESULTS AND DISCUSSION 



The viscous inviscid interaction code was run with several cascades on which ex- 
perimental data is available. In order to enable a thorough comparison between exper- 
imental results and the computed results, a very detailed experimental data base is 
needed. The data should include measurements of the boundary layer development along 
the blade, velocity profiles along the boundary layer, transition location and distribution, 
flow separation, and external velocity distribution. 

Unfortunately, very few cascade experiments has been performed, which obtained 
the required data with sufficient accuracy, due mostly to the lack of appropriate meas- 
urement equipment. Only recently, with the introduction of non--interfering methods 
like the Laser Doppler Velocimeter (LDV), the required data can be measured accu- 
rately. 

Recently an experiment involving the investigation of a linear compressor cascade 
of Controlled Diffusion Blading (which will be referred here as the CD cascade) has been 
carried out by Elazar [Ref. 8]. Most of the work in the present thesis, involves compar- 
ison of the computer code results with Elazar's experimental results. 

Other cascades that were investigated are: 

1. A shockless, supercritical airfoil cascade, designed in 1974 by Korn in cooperation 
with Pratt & Whitney Aircraft (referred here as the P & W cascade). The exper- 
imental results of the cascade were obtained from a report by Hobbs, Wagner, 
Dannenhoffer and Dring [Ref. 9]. 

2. Stator blade of a single stage axial compressor (referred here as the C4 cascade). 
The blade profile is the British C4 section (10% thickness) on a circular arc camber 
line. The experiment has been performed by Walker [Ref. 10}. The detailed 
boundary layer measurements are not presented in the report and were obtained 
directly from the author. 

The code failed to run with two other cascades: 

1. A highly loaded, double circular arc blade with a sharp leading edge and a sharp 
trailing edse, used in a compressor cascade that was investigated by Deutsch and 
Zierk [Ref. 11], 

2. V2 double circular arc blade, highly loaded cascade. This cascade was investigated 
by Hoheisel and Seyb [Ref. 12]. 

In both cases the code calculated the potential flow successfully but failed in trying 
to compute the first cycle of the boundary layer calculations. 
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A. CD CASCADE 



The experimental data for the CD cascade was obtained at M = 0.25. R, = 700000 
and at three inlet angles: 40° (the design condition). 43.4° and 46° . The spacing was 

0.6 of the chord and the stagger angle 14.27° . A general layout of the cascade F shown 
in Figure 1 on page 20. 

The following observations were concluded from the experiment: 

1 . A separation bubble exists near the leading edge on the upper surface at all the inlet 
angles. The bubble became larger at increased inlet angles. 

2. Transition from laminar to turbulent flow occurred above the separation bubble 
(on the upper surface). 

3. Transition on the lower surface occurred at midchord. 

4. The boundary layer thickness on the upper surface increased with inlet angle, and 
reached a thickness of 15% chord at the highest inlet angle. The boundary layer 
thickness on the lower surface did not chanee sienificantlv with inlet anele. 

5. The turbulent boundary layer on both surfaces remained fully attached at all the 
inlet angles. 

1. Transition location and intermittency distribution 

The effects of the transition location and the intermittency factor were investi- 
gated. The code was first run with the transition location calculated by the code, and 
with several values of the intermittency factor G y . It was found that the code did not 
run with G y = 1200 (which is the value used usually for high Reynolds numbers). The 
highest value of G y with which the code run successfully was 900. 

The code failed to predict the separation bubble on the upper surface, and pre- 
dicted laminar separation at 78% chord on the lower surface (which did not occur in the 
experiment). Transition on the upper surface occurred at 41% chord (detected by 
Michel's criterion) and at 78% chord on the lower surface (at laminar separation). 

The shape factor computed by the code was compared to the experimental re- 
sults. As can be seen in Figure 2 on page 21 the shape factor as predicted by the code 
deviates substantially from the actual results, due mainly to the different transition lo- 
cation. 

On the lower surface, as can be seen in Figure 3 on page 22 the shape factor 
deviates even more from the experimental results. In this figure the effect of changing 
the intermittency factor G y can be seen. For both the extreme values of G y . 10 and 900, 
the computed shape factor curve is far from agreement with the actual results. 
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Figure 1. Controlled Diffusion cascade 
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Figure 2. Shape factor comparison on the upper surface: Transition computed by 

the code (/? = 40°) 
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Figure 3. Shape factor comparison on the lo>>er surface: Transition computed by 

the code (/? = 40°). 
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Since the code did not predict the existence of the separation bubble on the 
upper surface, it was decided to try to eliminate the smoothing of the external velocity 
as computed by the potential flow subroutine (originally, the velocities were smoothed 
in subroutine CASBLP prior to boundary layer calculations). It was found that without 
smoothing the velocities a small separation bubble is predicted by the code at 4% of 
chord. The onset of transition is set by the code at the beginning of the separation 
bubble. 

The code was run with unsmoothed velocities with two values of G y , 10 and 900. 
The shape factor behavior can be seen in Figure 4 on page 24. Changing the value of 
G v did not change the shape of the curve much, and generally the shapes of the com- 
puted and the experimental curves look alike. 

The elimination of the velocitv smoothing in the code, also affects the thickness 
of the boundary layer. In Figure 5 on page 25 the displacement thickness is plotted for 
both cases (with and without the velocity smoothed). Without smoothing, the displace- 
ment thickness is much thicker, especially on the rear half of the blade, which is closer 
to the actual results. 

The effect of changing the intermittency distribution to the one used by Rodi 
and Schonung [Ref. 7] was investigated. It was found, as can be seen in Figure 6 on page 
26 that the effect of the new model is equivalent to using G y in the present model. 

On the lower surface it was necessary to run the code with transition as input, 
to get reasonable results, as can be seen in Figure 7 on page 27 for transition input at 
21% of chord. 

At the oir design conditions (inlet angles of 43.4° and 46°) a similar behavior of 
the transition has been observed, as can be seen for example in Figure 8 on page 28 for 
the upper surface and in Figure 9 on page 29 for the lower surface, both at /? = 46°. 
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Figure 4. Shape factor on the upper surface without velocity smoothing. 
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Figure 5. Displacement thickness: The effect of velocity smoothing. 
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Figure 6. The effect of the intermittency model: Upper surface, (i = 40° 
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Figure 7. Shape factor on the loner surface with transition input at 21%. 
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Figure 8. Shape factor at P = 46° on the upper surface 
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Figure 9. Shape factor at /J = 46° on the loner surface. 
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2. External Velocity Distribution 

The external velocity distribution, computed by the code using the interaction 
law, was compared to the experimentally measured velocities. It was found that in gen- 
eral, the velocities measured experimentally, were higher than those computed by the 
code for all the inlet angles. 

There are two possible sources to the discrepancy in the velocities: 

1. The computer code calculates pure 2— D flows. In the experiment the flow was 
observed to accelerate due to the effect of the boundary layer on the side walls (a 
3--D effect). This effect was calculated in the experiment and is referred to as the 
AVDR correction [Ref. 8, p.43]. 

2. The flow accelerates due to the thickening of the boundary layer. Since the 
boundary layer as computed by the code is substantially thinner than the actual 
boundary layer (as will be discussed in the next section) the external velocities 
predicted by the code are smaller. 

To compensate for the first error source, all the computed velocities were com- 
pared to the experimental velocities corrected by the AVDR correction. The comparison 
between the velocities can be seen in Figure 10 on page 31 for /? = 40°, in Figure 11 
on page 32 for /? = 43.4° and in Figure 12 on page 33 for /? = 46°. 

It can be seen from the figures that the difference between the computed and the 
experimental velocities is larger on the lower surface. The reason might be the method 
with which the correction to the inviscid velocity is computed. The assumption on which 
the interaction law is based, is that only sources (representing the viscous effects) on the 
surface being considered, affect the local velocity. In reality, the boundary layer on both 
surfaces affects the local velocity (because the boundary layer developed on the upper 
surface of a blade, causes a velocity disturbance that is felt on the lower surface of the 
adjacent blade). 

Since the boundary layer on the lower surface is much thinner, its effect on the 
velocity on the upper surface is much smaller than the effect of the upper surface 
boundary layer on the lower surface velocity. 
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Figure 10. External velocity at /? = 40° 



31 



o 




Figure 11. External velocity at /? = 43.4° 
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Figure 12. External velocity at fi = 46° 
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3. Boundary Lajer Thickness 

The boundary layer thickness as computed by the code was compared with the 
experimental results by comparing the displacement thicknesses. 

It was found that on the lower surface the computed and the actual displace- 
ment thickness agree quite well, as can be seen in Figure 13 on page 35 for /? = 40°, 
Figure 14 on page 36 for /? = 43.4° and in Figure 15 on page 37 for ft = 46° . 

On the upper surface the displacement thicknesses computed by the program 
are significantly thinner than those measured experimentally. The difference between the 
computed and the actual thickness increases along the blade and it increases with in- 
creased inlet angle. It was found that by using a different expression for the inner region 
eddy viscosity (as mentioned in chapter II), the displacement thickness can be increased, 
but the difference between the actual and the computed thickness is still substantial, es- 
pecially at the higher inlet angles. Figure 16 on page 38, Figure 17 on page 39 and 
Figure 18 on page 40 shows the displacement thickness on the upper surface for the 
three inlet aneles, with the original and the modified eddv viscosity models. 

The large error in the prediction of the boundary layer thickness, can be the re- 
sult of several reasons: 

1. The transition model used in the code, sets the onset of transition at the first point 
of laminar separation. It causes rapid transition to turbulent flow which reattaches 
immediately, resulting in a very small separation bubble compared to the bubble 
observed in the experiment. 

2. The turbulent model used in the code could be inaccurate. It was derived based on 
empirical data obtained in single airfoil experiments and not with cascades. In ad- 
dition the present model does not include the effects of the free stream turbulence 
(that was relatively high in the experiment). 

3. The boundary layer as measured in the experiment was quite thick, especially at the 
higher inlet angles (it reached 15% of the chord at /? = 46°). Such a thick bound- 
ary layer may violate the basic assumptions on which the boundary layer 
equations, and the interaction law, were based (especially when the spacing be- 
tween the blades is small, 60% chord in this case). 

It was suggested that one of the possible reasons to the inaccurate prediction 
of the boundary layer is the blunt trailing edge of the blade, that might cause difficulties 
in the computations. A modified blade, with a sharp trailing edge has been run, and the 
displacement thickness distribution can be seen in Figure 19 on page 41. As can be seen 
in the figure the sharp trailing edge affects only the boundary layer adjacent to the 
trailing edge, and therefore cannot provide an explanation to the difference between the 
actual and the computed displacement thickness. 
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Figure 13. Displacement thickness on the loner surface (/l = 40° ) 
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Figure 14. Displacement thickness on the loner surface ((1 = 43.4° ) 
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Figure 15. Displacement thickness on the louer surface (// = 46° ) 
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Figure 16. Displacement thickness on the upper surface (ft = 40° ) 
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Figure 17. Displacement thickness on the upper surface (fl = 43.4° ) 
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Figure 18 . Displacement thickness on the upper surface (/J = 46 °) 
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Figure 19. The effect of sharp trailing edge. 
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4. Comparison (o a Navier Stokes Code 

A limited comparison of the experimental and the computed results with a 
Navier Stokes (N.S.) code calculations has been performed. The N.S. code has been de- 
veloped and run by S. J. Snamroth of Scientific Research Associates inc. in cooperation 
with Pratt and Whitney Aircraft. 

Since the N.S. code does not compute the displacement thickness, the velocity 
profiles near the surface of the blade were compared. The comparisons were made at 
90% chord on the suction surface for all three inlet angles. 

At the design point. /? = 40°, shown in Figure 20 on page 43. both the inter- 
active code and the N.S. code failed to predict accurately the actual velocity profile. In 
this case the interactive code seems to yield somewhat better results than the N.S. code. 

At the higher inlet angles. /? = 43.4° and /? = 46°, shown in Figure 21 on page 
44 and in Figure 22 on page 45 respectively, the N.S. calculations show significantly 
better agreement with the experimental results than the interactive code. 

From these comparisons, it can be seen that the interactive code deviation from 
the actual results increases with increased inlet angle (increased loading of the cascade), 
whereas the N.S. code deviation seems to decrease with increased inlet angle. 
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Figure 20. Die results of the N. S. code at p = 40° 
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Figure 21. The results of the N. S. code at /? = 43.4° 
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Figure 22. The results of the N. S. code at (i = 46° 



B. P & W CASCADE 



The experimental data for the P &. W cascade was obtained at inlet flow angle of 
52°, at M = 0. 11, and Reynolds number of 478000. The cascade had a stagger angle of 
15.75° and 0.7 spacing. A general layout of the cascade is shown in Figure 23 on page 
47. 

A comparison of the computed and the measured pressure coefficients on the blade 
is shown in Figure 24 on page 48. There is a good agreement between the computed 
and the measured C P . 

The displacement thickness was measured in the experiment only at 96.8% of chord. 
This measurement is compared to the computed results in Figure 25 on page 49. As can 
be seen, the computed and the measured data agree almost perfectly on the lower sur- 
face, and quite well on the upper surface. The difference observed on the upper surface 
is caused by the early prediction, by the code, of trailing edge separation, a short dis- 
tance upstream of the actual location. This can also be observed when comparing the 
velocity profiles at that point, in Figure 26 on page 50. The computed velocity curve 
shows a small zone of reversed flow near the surface of the blade. This reversed flow 
could be the result of a too early prediction of trailing edge separation by the code, or 
it could have existed in the actual flow but not detected because of its size. 
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P&W CASCADE 





Figure 23. Pratt & Whiteny cascade 
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Figure 24. Comparison of pressure coefficient. 
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Figure 25. Displacement thickness comparison: Experimental data shown at 

96.8% chord. 
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Figure 26. Velocity profile at 96.8% chord on the upper surface. 
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C. C4 CASCADE 



The C4 cascade has a stagger angle of 29.5°, a camber angle of 31.1° and spacing 
of 0.992. It has been tested at Reynolds numbers of about 200000. and inlet angles of 
34.1° to -17.7° (which corresponds to incidence angles of -10.9° to 2.7 C ). The general 
layout of the cascade is shown in Figure 27 on page 52. A computer code that generates 
the coordinates of the blade and a summary of some experimental results arc given in 
Appendix B. 

The code was run with the intcrmittcncy constant G y = 10. Higher values of G v 
(abose 100) caused numerical problems in the code. The onset of transition was first 
taken at the point where it was observed in the experiment. At the lower inlet angle it 
seems that a better agreement with the experimental results can be obtained by delaying 
the onset of transition but trying to implement it resulted in numerical breakdown of the 
computation. At the higher inlet angles, better agreement with the experimental results 
was achieved by initiating the transition earlier (at 26% chord for p = 45.6° and at 21% 
for /? =47.7° as compared to 44% and 36% chord as observed in the experiment). 

1. Displacement Thickness 

Comparisons of the experimental data to the computed displacement thickness 
are shown in Figure 2S on page 53 for inlet angle of 34.1°. in Figure 29 on page 54 for 
inlet angle of 36.3°, in Figure 30 on page 55 for inlet angle of 45.6° and in Figure 31 
on page 56 for inlet angle of 47.7°. 

As can be seen in the figures, there is a good agreement between the actual and 
the computed results at the two lower angles (p = 34.1° and/? = 36.3°, in which the 
incidence angles were negative). At the two higher angles, p = 45.6° and p = 47.7° the 
computed results agree with the actual results up to about 70% chord, and then the 
displacement thickness predicted by the code becomes much thicker than the actual one. 

The code predicted a large flow separation area starting at about 70% chord at 
the lower inlet angles, and at about 46% chord at the higher inlet angles. This flow 
separation was not obsc-ved in the experiment. The discrepancies between the computed 
and the actual results behind 60% to 70% chord can be explained by the inaccurate 
calculations by the code due to the large separated areas. When the code encounters 
separation, several approximations are made (like the FLARE approximation) based on 
the assumption that the separated area is small. When the separated area is large, these 
approximations may result in inaccurate prediction of the flow field. 
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Figure 27. C4 Cascade 
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Figure 28. C4 cascade at /? = 34.1°: Displacement thickness comparison with 

computed results. 
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Figure 29. C4 cascade at ft = 36.3°: Displacement thickness comparison with 

computed results. 
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Figure 30. C4 cascade at [i = 45.6°: Displacement thickness comparison with 

computed results (G= 10). 




Figure 31. C4 cascade at /? = 47.7°: Displacement thickness comparison with 

computed results (G= 10). 
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2. External Velocity and Velocity Profiles Comparisons 

A comparison of the external velocity on the upper surface of the blade is shown 
in Figure 32 on page 5S for inlet angle of 45.6° and in Figure 33 on page 59 for inlet 
angle of 47.7°. h can be seen that there is a good agreement between the experimental 
and the computed results up to about S0% chord. Near the trailing edge the computed 
results deviate from the experimental results due to the inaccuracy in the calculations 
of the displacement thickness. 

A comparison of the velocity profiles in the boundary layer at 50% chord is 
shown in Figure 34 on page 60 for inlet angle of 34.1° and in Figure 35 on page 61 for 
inlet angle of 36.3°. The agreement between the calculated velocity profiles and the 
measured \elocity profiles is very good. 
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Figure 32. C4 cascade at P = 45.6°: L:\ternal velocity distribution. 
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Figure 33. 



C4 cascade at / 1 - 47 . 7 °: External velocity distribution. 



CO 




Figure 34. C4 cascade at /J = 34.1°: Velocity profile at 50% chord. 
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Figure 35. C-4 cascade at ft = 36.3°: Velocity profile at 50% chord. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The interactive viscous inviscid computer code, has been investigated by comparing 
its predictions of boundary layer parameters to experimental data. 

It has been found that the code yields reasonable results for lightly loaded cascades, 
but the prediction of the boundary layer thickness on the suction surface of highly 
loaded cascades deviates significantly from experimentally measured data. In two cases 
involving highly cambered cascade blades, with sharp leading edge, the code failed to 
run. It was also found that the prediction of external velocity distribution on highly 
loaded cascades was inaccurate. 

The main reasons to the discrepancies in the prediction of the boundary layer 
thickness seem to be: 

1. Inaccuracy in predicting How parameters in regions of large flow separation (due 
to inadequate transition model and approximations made in calculating the How in 
separated areas). 

2. Inaccurate turbulence modelling. 

3. Possible violation of the basic assumptions of the boundary layer theory in areas 
of very thick boundary layers. 

4. The wake is not calculated by the code. The result is inaccurate flow prediction 
near the trailing edge. 

The inaccuracy in the prediction of the external velocity distribution in highly loaded 
cascades is due to the interaction law. which docs not account for the presence of adja- 
cent blades. 



B. RECOMMENDATIONS 

The recommended steps in order to improve the code are: 

1. Improving the interaction law by assuming a distribution of sources on the actual 
surface (instead of the assumption of a flat plate), letting the correction term to the 
external velocity vary across the boundary layer and distributing sources on the 
adjacent blade as well for better modelling of the boundary layer effect on the ex- 
ternal velocity. 
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2. Changes to the derivation of the boundary layer equations should be investigated 
to allow a better treatment of thick boundary layers t like omitting the assumption 
of cpidy — 0 across the boundary layer). 

3. Different turbulence models should be investigated. 

4. The wake should be included in the calculations. 
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APPENDIX A. COMPUTER CODE LISTING 



C*** 

c*** 
c*** 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
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c 
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c 
c 
c 
c 
c 
c 
c 
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VISCOUS-INVISCID INTERACTION PROGRAM FOR CASCADE FLOWS 



.-***** I NT 0 0 0 1 
***** INTO 002 



3. A 



*********** ”*"**<?********* <c** INTO 003 

INT0004 
INT0005 
INT0006 
INT0007 

THIS VISCOUS-INVISCID INTERACTION METHOD, CAPABLE OF COMPUTING BOTH INT0008 



VERSION 
JANUARY 87 



SINGLE AIRFOIL AND CASCADE FLOWS, WAS DEVELOPED BY CEBECI AND 
COLLABORATEURS AT LONG BEACH STATE AND DOUGLAS AIRCRAFT COMPANY. 
THE CODE APPLIES TO INCOMPRESSIBLE, 2-DIMENSIONAL, STEADY FLOWS 
PAST LINEAR, ARBITRARILY STAGGERED CASCADES. THE METHODS BASIC 
INGREDIENTS INCLUDE 

1. A FIRST ORDL1 PANEL METHOD TO SOLVE LAPLACE'S EQUATION, 



INT0009 
INT0010 
INT0011 
INT0012 
INTOO 13 
INT0014 



2. A FINITE DIFFERENCE SCHEME TO SOLVE THE BOUNDARY LAYER EQUATIONS INTOO 15 



3. 



SUBJECT TO DIRECT OR INTERACTIVE BOUNDARY CONDITIONS, 

A STRONG INTERACTION MODEL TO COUPLE VISCOUS AND INVISCID FLOW 
RESULTS, AND 

A ZERO EQUATION, ALGEBRAIC TURBULENCE MODEL TO ESTIMATE 
TURBULENT SHEAR STRESSES. 



INT0016 
INT0017 
INT0018 
INTOO 19 
INT0020 
INT0021 

IN SUMMARY, THE CODE WILL PROVIDE, FOR ATTACHED AS WELL AS MODERATE -INT0022 
LY SEPARATED FLOWS PAST SINGLE AIRFOILS OR CASCADES, THE FOLLOWING INT0023 
1. INVISCID AND VISCOUS PRESSURE DISTRIBUTIONS, INT0024 

DISTRIBUTIONS OF INT0025 1 

A. LOCAL SKIN FRICTION COEFFICIENT, INT0026' 

B. DISPLACEMENT AND MOMENTUM THICKNESS, AND * INT0027> 

VELOCITY PROFILES ACROSS THE BOUNDARY LAYER. INT0028' 

INT0029 1 

MODIFICATIONS SINCE VERSION 3. 0: INT0030i 

1. PRECISE ASSIGNMENT OF BEGIN OF TRANSITION. INT0031' 

2. CORRECTION OF AN ERROR IN THE CALCULATION OF MOMENTUM THICKNESS. INT0032' 

3. ADDITIONAL PRINT OPTION: IP=-2 WILL PROVIDE AN INPUT FILE (UNIT INT0033' 



2 . 



3. 



NUMBER 12) FOR THE PLOTTING ROUTINE. 



COMMON /BLCO/ NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 
COMMON/BLOW/ VN( 100) 

COMMON/BLIN/ TITLE(20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR , ISTRP , I CYCLE , ICYTL, XCTRS( 2) ,TRFIND( 2 ) 

COMMON/C ASCDE/ INLET , SP , SINGLE , ALPHAA , ALPHAI , STAG 
COMMON/TRN/ PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON/PLOT/NVPC 2) ,NXVP( 20 , 2 ) , ICC 

DIMENSION XO( 100) , YO( 100) ,X( 100) ,Y( 100) , VCOM( 100) ,DLS( 100) , 
+ XS( 100) , YS( 100) , XSTGR( 100) , YSTGR( 100) ,DBPP( 100) 

DIMENSION CASEID( 20) ,XCTRI( 2) , ITRI( 2) ,NBL(2) 

LOGICAL SINGLE, TRFIND 
TRFIND( 1)= .FALSE. 

TRFIND( 2 )= .FALSE. 



INT0034 1 

INT0035 1 

INT0036 1 

INT0037I 

INT0038 1 

INT0039 1 

INT0040 1 

INT0041 1 

INT0042 1 
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INT0044 
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J.CASE = 0 


INT00500 


1 


READ(5,5,END=999) TITLE 


INT005 10 


5 


FORM AT ( 20A4) 


INT00520 




ICASE = ICASE + 1 


INT00530 




REWIND 3 


INT00540 


C 


READ (5,10) 


INT00550 


10 


FORM AT ( IX) 


INT00560 




READ (5,20) ITRI(1),ITRI(2), IRST, ICYTL, IP 


INT00570 


20 


FORMAT( 1615) 


INT00580 


C 


READ (5,10) 


INTO 05 90 




READ (5,25) INLET, I STAG, ALPHAI , STAG, SP,PGAMTR, OMEGA 


INT00600 


25 


FORMAT( 215 ,5F10. 0) 


INT006 10 




READ (5 ,27)RN,XCTRI( 1) ,XCTRI(2) , ALPHA A 


INT00620 


27 


FORMAT( 4E10. 0) 


INT00630 




IF (IP. EQ. -2) THEN 


INT00640 




READ (5,20) NVP( 1) ,NVP( 2) 


INT00650 




IF (NVP(l).NE. 0) READ (5,20) (NXVP( I , 1) , 1=1 ,NVP( 1) ) 


INT00660 




IF (NVP(2). NE. 0) READ (5,20) (NXVP( 1,2), 1=1 ,NVP( 2) ) 


INT00670 




END IF 


INT00680 




IF (ICASE .EQ. 1) READ (5,20) N,NI 


INT00690 




I READ = 1 


INT00700 




I BLOW = 1 


INT007 10 




SINGLE = .FALSE. 


INT00720 




IF (SP . LE.- 0.0) SINGLE = .TRUE. 


INT00730 




N = N - 1 


INTO 07 40 




Nl= N + 1 


INT00750 




IF (ICASE .GT. 1) THEN 


INT00760 




N1 = N1SAVE 


INT007 7 0 




N = N - 1 


INT00780 




GOTO 53 


INT00790 




END IF 


INT00800 




IF (IREAD .EQ. 1) GO TO 40 


INT00810 


C 


READ (5,10) 


INT00820 




READ ( 5 , 3C ) ( XO( I ) , YO(I), 1=1, N+l) 


INT00830 


30 


FORMAT( 2F10. 0) 


INTO 0840 




GO TO 50 


INT00850 


C 




INT00860 


C40 


READ( 5,10) 


INT00870 


40 


READ(5 ,45) (XO(I) , 1=1, N+l) 


INT00880 


C 


READ( 5,10) 


INT00890 




RE AD (5 ,45) (YO(I) , 1=1, N+l) 


INT00900 


45 


FORM AT (6F 10. 0) 


INT00910 


C 




INT00920 


50 


CONTINUE 


INT00930 




IF (IP. EQ. -2) THEN 


INT00940 




WRITE( 12,20) N+1,NVP( 1) ,NVP(2) ,9 0,70, INLET 


INT00950 




IF (NVP(l).NE. 0) WRITE( 12,20) (NXVP(I , 1) ,I=1,NVP( 1)) 


INT00960 




IF (NVP( 2) . NE. 0) WRITE( 12,20) (NXVP( I , 2) , 1=1 ,NVP( 2) ) 


INT00970 




IF ( INLET. NE. 1) WRITE( 12,80) RN , ALPHAA 


INT00980 




IF ( INLET. EQ. 1) WRITE(12,80) RN, ALPHAI 


INT00990 




WRITE( 12,82) (X0( I ) , 1=1 , N+l) 


INTO 1000 




WRITE( 12,82) ( Y0( I ) , 1=1 ,N+1) 


INT01010 


80 


F0RMAT(2E15. 5) 


INTO 1020 


82 


FORMAT( 8F10. 6) 


INTO 1030 




END IF 


INTO 1040 


C 




INTO 1050 
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NRITE = (NI+l)/2 

IMIN = (N 1-1)/ 2+1 

IF( (Nl/2*2) . EQ. Nl) IMIN = Nl/2 

CALL TRGRID (N1,X0, Y0,NI , NRITE, 0. 5 , IMIN, RAD, 1 ,NXSS1) 

N1SAVE = Nl 
53 CONTINUE 

ALPHAA = 0.0174533 * ALPHAA 
ALPHAI = 0.0174533 * ALPHAI 
STAG = 0.0174533 * STAG 
C 

IF (INLET .EQ. 0) THEN 

ALPHA = ALPHAA 

ELSE 

ALPHA = ALPHAI 
END IF 
C 

IF (ISTAG .NE. 0) THEN 

CALL STAGR( N 1 , STAG , XO , YO , XSTGR , YSTGR ) 

ELSE 

DO 55 I = 1 , Nl 
XSTGR(I) = XO(I) 

YSTGR(I) = YO( I ) 

55 CONTINUE 

END IF 

READ DATA FROM VISCOUS CAL. 

ICYCLE = 0 
60 ICYCLE = ICYCLE + 1 

C 

CALL POTNL( N 1 , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , VCOM , 

+ DBPP) 

IF (ICYCLE . GT. ICYTL) THEN 
REWIND 3 

WRITE ( 3)N1 , (XO( I ) , YO( I ) ,DLS( I ) , VN( I) ,DBPP( I ) , 1=1 ,N1) 

GOTO 1 
END IF 
C 

IF (ISTAG .NE. 0) THEN 

DO 70 I = 1 , Nl-1 

X(I) = 0.5 * (XO(I)+XO(I+l)) 

Y( I ) = 0.5 * ( YO( I )+YO( 1+1) ) 

70 CONTINUE 
END IF 
C 

CALL CASBLP(Nl,XO,YO,X,Y,XS,YS, DLS, VCOM, DBPP, RN 
+ ,NBL, ITRI ,XCTRI , TITLE) 

GO TO 60 
999 CONTINUE 
STOP 
END 
C 

SUBROUTINE POTNL( Nl , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , 
+ VCOM, DBPP) 

C 

COMMON/BLOW/VN( 100) 
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INTO1350 
INTO 13 60 
INTO1370 
INTO 1380 
INTO1390 
INT01400 
INTO 1410 
INTO1420 
INTO1430 
INTO 1440 
INTO 1450 
INTO1460 
INTO 1470 
INTO 1480 
INTO1490 
INT01500 
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INTO1520 
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INTO1550 
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INTO 15 90 
INT01600 
INTO 16 10 
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COMMON /BLCO/ NX ,NXT ,NP ,NPT,NTR, IT, INVRS ,NS , IP 

COMMON/ BLIN/TITLE( 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR , I STRP , IC YCLE , IC YTL , XCTRS ( 2 ) , TRF IND ( 2 ) 

COMMON/ CASCDE / INLET, SP , SINGLE , ALPHAA , ALPHA I , STAG 
C SIMPLE SOURCE POTENTIAL CODE 

DIMENSION AOFF( 100, 100) , BOFF( 100 , 100) , XP(IOO), YP( 100) ,X( 100) , 
+ S(100) ,C(100), D( 100) , VTAN( 3,100) , VNOR( 3, 100), R( 3, 100), 

+ VCOM( 100) , SIGCOM( 100) ,CP( 100) ,X0( 100) ,Y0( 100) 

+ , VNC( 100) ,D1( 100) ,D2( 100) ,D3( 100) ,S0( 100) ,SC( 100) 

+ ,XOFF( 100) , YOFF( 100) ,T(3 , 100) , VTCOM( 100) , VNCOM( 100) 

+ ,XS( 100) ,YS( 100) ,SOFF( 100) ,COFF( 100) ,XPOFF( 100) ,YPOFF( 100) 

+ ,Y( 100) ,SIG(3, 100) ,DLS( 100) ,DLSC( 100) ,A(100, 100) ,B( 100,100) 

+ , XSTGR( 100) ,YSTGR( 100) , VUT(3) ,VLT(3) , VUN(3) ,VLN(3) ,DBPP( 100) 

+ ,CPI( 100) ,XOS( 100) , YOS( 100) ,DBPPC( 100) 

REAL NUM1 , NUM2 
LOGICAL OFF, SINGLE 
OFF = .FALSE. 

PI = 3. 1415S2 
CM = 0. 0 
N = N1 - 1 

IF (ICYCLE . EQ. 1) THEN 
IF (IRST .EQ. 0) THEN 
DO 10 1=1, N1 
DLS(I) = 0. 0 
VN (I) = 0. 0 
DBPP( I )= 0. 0 
10 CONTINUE 

ELSE 

DO 5 I = 1 , N1 

XOS(I) = XO(I) 

YOS(I) = Y0( I) 

5 CONTINUE 

READ (3) NT, (XS( I) , YS( I) ,DLSC( I) , VNC( I) ,DBPPC( I) , 1=1 , NT) 

XMIN = XS(1) 

DO 15 I = 2 , NT 

IF (XS(I) .GT. XMIN) GOTO 15 

XMIN = XS(I) 

IMIN = I 

15 CONTINUE 

DO 17 I = 1 , NT 
IF (I .LT. IMIN) GOTO 16 
XS(I) = XSCI) - XMIN 
GOTO 17 

16 XS(I) = XMIN - XS(I) 

17 CONTINUE 
C 

XMIN = XOS(I) 

DO 20 I = 2 , N1 

IF (XOS(I) .GT. XMIN) GOTO 20 

XMIN = XOS(I) 

IMIN = I 
20 CONTINUE 

DO 22 I = 1 , N1 

IF (I .LT. IMIN) GOTO 21 

XOS(I) = XOS(I) - XMIN 



INT01620 
INTO 1630 
INTO 1640 
INTO 1650 
INTO 1660 
INTO 16 70 
INTO 1680 
INT01690 
INT01700 
INT017 10 
INT01720 
INT01730 
INT01740 
INTO 1750 
INT01760 
INT01770 
INTO 17 80 
INTO 1790 
INTO 1800 
INT01810 
INTO 1820 
INTO 1830 
INT01840 
INTO 1850 
INTO 1860 
INT01870 
INTO 1880 
INT01890 
INTO 19 00 
INTO 19 10 
INT01920 
INT01930 
INT01940 
INT01950 
INTO 19 60 
INT01970 
INTO 1980 
INT01990 
INT02000 
INT02010 
INT02020 
INT02030 
INT02040 
INT02050 
INT02060 
INT02070 
INT02080 
INT02090 
INT02100 
INT02110 
INT02120 
INT02130 
INT02140 
INT02150 
INT02160 



67 



O ho 10 



GOTO 22 

1 XOS(I) = XMIN - XOS(I) 

2 CONTINUE 

CALL D IFF3 ( NT , XS , DLSC , D 1 , D2 , D3 , 0 ) 

CALL INTRP3 ( NT , XS , DLSC , D 1 , D2 , D3 , N 1 , XOS , DLS ) 

CALL AMEAN ( 1 ,N1 ,XOS ,DLS, 1) 

CALL DIFF3 (NT ,XS , VNC ,D1 ,D2 ,D3 , 0) 

CALL INTRP3 ( NT , XS , VNC , D 1 , D2 , D3 , N 1 , XOS , VN) 

CALL AMEAN ( 1 ,N1 ,XOS , VN, 1) 

CALL DIFF3 ( NT.XS , DBPPC ,D1 ,D2 ,D3 , 0) 

CALL INTRP3 ( NT , XS , DBPPC , D 1 , D2 , D3 , N1 , XOS , DBPP ) 

CALL AMEAN ( 1,N1, XOS, DBPP, 1) 

END IF 
END IF 
DO 30 1=1 ,N1 
XP(I) = XSTGR(I) 

YP( I ) = YSTGR(T) 

30 CONTINUE 

C CALCULATE GEOMETRIC QUANTITIES 
DO 100 J=1,N 

VNC(J) = 0.5 * (VN(J) + VN( J+l) ) 

X(J)= . 5*(XP(J)+XP( J+l)) 

Y( J)= . 5*(YP(J)+YP(J+1)) 

D(J)= SQRT((XP(J+1)-XP(J))**2 + ( YP( J+l) -YP( J) )**2) 

C(J)= (XP( J+l) -XP( J) )/D( J) 

S(J)= ( YP( J+l) -YP( J) )/D( J) 

100 CONTINUE 
C 

IF ( INLET .NE. 0 .AND. .NOT. SINGLE) THEN 
SUM = D(l) 

DO 35 J = 2, N 
SUM = SUM + D( J) 

35 CONTINUE 

Q = 2. 0 * PI * SUM / SP 
ELSE 
Q = 0. 0 
END IF 
C 

C CALCULATE NORMAL AND TANGENTIAL MATRICES 
102 CONTINUE 

IF (SINGLE) THEN 

IF ( .NOT. OFF) THEN 

DO 120 1=1, N 

DO 110 J=1,N 

IF (J . EQ. I) GO TO 105 

XX= (X(I)-X(J))*C(J) + (Y(I)-Y(J))*S(J) 

YY=-(X( I ) -X( J) )*S( J) + (Y(I)-Y(J) )*C( J) 

UU= LOG( ( (XX+. 5*D( J))**2+YY**2)/((XX-. 5*D( J) )**2+YY**2) ) 
VV= 2. *ATAN2( YY*D( J) , XX** 2+YY**2-( . 5*D( J) )**2) 

SS= S(I)*C(J) - C(I)*S(J) 

CC= C(I)*C(J) + S(I)*S(J) 

A( I , J)= -UU*SS + VV*CC 
B( I , J)= UU*CC + VV*SS 
GO TO 110 

105 A( I , J) = 6.2831853 
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B( I , J) = 0. 0 


INT02730 


no 


CONTINUE 


INT02740 


120 


CONTINUE 


INTO 2 750 




ELSE 


INT02760 




DO 140 1=1, N 


INT02770 




DO 130 J=1,N 


INT02780 




XX= (XOFF( I ) -X( J) )*C( J) + ( YOFF( I)-Y(J))*S(J) 


INT02790 




YY=-(XOFF(I)-X(J))*S(J) + (YOFF(I)-Y(J))*C(J) 


INT02800 




UU= LOG( ( (XX+. 5*D( J) ) ** 2+YY** 2 ) / ( ( XX - . 5*D( J) )**2+YY**2) ) 


INT028 10 




VV= 2. *ATAN2( YY*D( J) , XX**2+YY**2-( . 5*D( J) )**2) 


INTO 2820 




SS= SOFF( I)*C( J) - COFF( I )*S( J) 


INT02830 




CC= COFF( I)*C( J) + SOFF( I)*S( J) 


INT02840 




AOFF( I , J)= -UU*SS + VV*CC 


INTO 2850 




BOFF( I , J)= UU*CC + W*SS 


INT02860 


130 


CONTINUE 


INT02870 


140 


CONTINUE 


INTO 2880 




END IF 


INTO 2 890 




ELSE 


INT02900 




IF ( .NOT. OFF) THEN 


INT029 10 




DO 50 1=1, N 


INTO 29 20 




DO 40 J=1,N 


INTO 29 30 




IF (J . EQ. I) GO TO 45 


INT02940 




XX= (X(I)-X(J))*C(J) + ( Y( I) -Y( J) )*S( J) 


INTO 29 50 




YY=-(X( I) -X( J) )*S( J) + ( Y( I) -Y( J) )*C( J) 


INTO 29 60 




DX1 = PI *(X( T )-XP( J)) / SP 


INTO 29 70 




DY1 = PI *(Y(: )-YP(J)) / SP 


INT02980 




DX2 = PI *(X(I)-XP(J+1)) / SP 


INT02990 




DY2 = PI *( Y ( I ) -YP( J+l) ) / SP 


INT03000 




R1SQ = (C0SNfPXi))**2 - ( COS( DY1) )**2 


INTO 30 10 




R2S0 = ( COSH( DX2) )**2 - ( COS( DY2) )**2 


INTO 30 20 




UU = LOG( R1SQ/R2SQ) 


INTO 30 30 




NUM1 = DX1 * COSH(DXl) * SIN(DYI) - 


INT03040 




+ DY1 * SINH(DXl) * COS(DYl) 


INT03050 




DNUM1= DX1 * SINH(DXl) * COS(DYl) + 


INT03060 




+ DY1 * COSH(DXl) * SIN(DYl) 


INTO 30 70 




NUM2 = DX2 * C0SH(DX2) * SIN(DY2) - 


INT03080 




+ DY2 * SINH(DX2) * C0S(DY2) 


INT03090 




DNUM2= DX2 * SINH(DX2) * C0S(DY2) + 


INT03100 




+ DY2 * COSH( DX2) * SIN(DY2) 


INTO 3 110 




EXV = 2.0 * AT AN 2 ( NUM2 , DNUM2 ) - 2.0 * ATAN2(NUM1 ,DNUM1) 


INT03120 




VV= 2. *ATAN2( YY*D( J) , XX**2+YY**2-( . 5*D( J) )**2) 


INTO 3 130 




VV = VV + EXV 


INTO 3 140 




SS= S(I)*C(J) - C(I)*S(J) 


INT03150 




CC= C(I)*C(J) + S(I)*S(J) 


INTO 3 160 




A( I , J)= -UU*SS + VV*CC 


INTO 3 170 




B( I , J)= UU*CC + VV*SS 


INTO 3 180 




GO TO 40 


INT03190 


45 


A( I , J) = 6. 2831853 


INT03200 




B( I , J) = 0. 0 


INT03210 


40 


CONTINUE 


INT03220 


50 


CONTINUE 


INT03230 




ELSE 


INTO 3 240 




DO 70 1=1, N 


INT03250 




DO 60 J=1 ,N 


INT03260 




XX= (XOFF( I) -X( J) )*C( J) + (YOFF(I)-Y(J))*S(J) 


INTO 3270 




YY=-(XOFF(I)-X(J))*S( J) + (YOFF(I)-Y(J))*C(J) 


INT03280 
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DX1 = PI *(XOFF( I) -XP( J) ) / SP 
DY1 = PI *(YOFF(I)-YP(J)) / SP 
DX2 = PI *(XOFF( I) -XP( J+l) ) / SP 
DY2 = PI *( YOFF( I ) -YP( J+l ) ) / SP 
R1SQ = ( COSH( DX1 ) )**2 - (COS(DYl) )**2 
R2SQ = (COSH(DX2) )**2 - (COS(DY2) )**2 
UU = L0G(R1SQ/R2SQ) 

NUM1 = DXl * COSH(DXl) * SIN(DYl) - 

+ DY1 * SINH(DXl) * COS(DYl) 

DNUM1= DXl * SINH(DXl) * COS(DYl) + 

+ DY1 * COSH(DXl) * SIN(DYl) 

NUM2 = DX2 * COSH(DX2) * SIN(DY2) - 

+ DY2 * SINH( DX2 ) * COS(DY2) 

DNUM2= DX2 * SINH(DX2) * COS(DY2) + 

+ DY2 * COSH(DX2) * SIN(DY2) 

EXV =2.0 * ATAN2(NUM2 ,DNUM2) - 2.0 * ATAN2(NUM1 ,DNUM1) 
VV= 2. * AT AN 2 ( Y Y V -'D ( J ) , XX**2+YY**2-( . 5*D( J) )**2) 

VV = VV + EXV 

SS= SOFF( I)*C( J) - COFF( I)*S( J) 

CC= COFF( J.)' V C( J) + SOFF( I)*S( J) 

AOFF( I , J)= -UU*SS + VV'-'CC 
BOFF( I , J)= UU*CC + VV*SS 
60 CONTINUE 

70 CONTINUE 

END IF 
END IF 

NORMAL AND TANGENTIAL COMPONENTS OF FUNDAMENTAL SOLUTIONS 

DO 160 1=1, N 
SUMR= 0. 

SUMT= 0. 

IF ( .NOT. OFF) THEN 
R( 1 , I )= S ( I ) +VNC ( I ) / COS ( ALPHA) 

T( 1 , I )= C( I) 

R( 2 , 1 )= -CCD 
T( 2 , 1 )= SCD 
DO 145 J=1,N 
SUMR = SUMR + B( I , J) 

SUMT = SUMT + A( I , J) 

145 CONTINUE 
ELSE 

R(1,I) = SOFF(I) 

T( 1 , 1 ) = COFF(I) 

R( 2 , 1 ) =-COFF( I ) 

T( 2 , I ) = SOFF(I) 

DO 150 J=1,N 

SUMR= SUMR + BOFF( I , J) 

SUMT= SUMT + AOFF( I , J) 

150 CONTINUE 
END IF 

R( 3 , I )= SUMR 
T ( 3 , I ) = SUMT 
160 CONTINUE 
C 

IF ( OFF ) GO TO 275 
C DECOMPOSITION OF MATRIX A 





DO 230 1=1, N-l 


INTO 3 850 




DO 220 K=I+1,\ T 


INT03860 




A(K,I)= A(K, I )/A( I j I) 


INT03870 




DO 210 J=I+1,N 


INT03880 




A(K, J)= A ( K , J ) - A(K,I)*A(I,J) 


INT03890 


210 


CONTINUE 


INT03900 


220 


CONTINUE 


INT039 10 


230 


CONTINUE 


INT03920 


C 1 


OPERATE ON FUNDAMENTAL RIGHT SIDES WITH LOWER TRIANGULAR 


INT03930 




DO 270 K=1 j 3 


INT03940 




DO 260 J=1 jN-1 


INT03950 




DO 250 I=J+1,N 


INT03960 




R(K j I)= R(K , I ) - A(I,J)*R(K,J) 


INT039 7 0 


250 


CONTINUE 


INT03980 


260 


CONTINUE 


INT03990 


270 


CONTINUE 


INT04000 


c : 


BACK SOLUTION 


INT04010 




DO 300 K=1 , 3 


INT04020 




DO 290 I=N,1,-1 


INT04030 




SUM= 0. 


INT04040 




DO 280 J=N, I-fl , -1 


INT04050 




SUM= SUM + A(I,J)*SIG(K,J) 


INT04060 


280 


CONTINUE 


INT04070 




SIG(K, I)= (R(K,I)-SUM)/A(I,I) 


INT04080 


290 


CONTINUE 


INT04090 


300 


CONTINUE 


INT04100 




OFF = .TRUE. 


INT04110 


C 




INT04120 


C 


ADD DIS -PLACE VERTICALLY TO THE BODY TO GENERATE 


INT04130 


C 


DISPLACEMENT SURFACE 


INT04140 


C 




INT04150 




DO 305 I = 2 , N 


INT04160 




SOFF(I) = (S(I)*D(I-1)+S(I-1)*D(I))/(D(I-1)+D(I)) 


INT04170 




COFF(I) = (C(I)*D(I-1)+C(I-1)*D(I))/(D(I-1)+D(I)) 


INT04180 


305 


CONTINUE 


INT04190 




SOFF(l) = 2. 0*S( 1) - SOFF( 2) 


INT04200 




S0FF(N1)= 2. 0*S(N) - SOFF(N) 


INT042 10 




COFF(l) = 2. 0*C(1) - C0FF(2) 


INT04220 




C0FF(N1)= 2. 0*C(N) - COFF(N) 


INT04230 




DO 306 I = 1 , N1 


INT04240 




XPOFF(I) = XP(I) - SOFF(I) * DLS(I) 


INT04250 




YPOFF(I) = YP( I ) + COFF(I) * DLS(I) 


INT04260 


306 


CONTINUE 


INT04270 




DO 307 I = 1 , N 


INT04280 




XOFF(I) = 0.5 * (XPOFF(I) + XP0FF(I+1)) 


INT0429 0 




YOFF(I) =0.5 * (YPOFF(I) + YP0FF(I+1)) 


INT04300 




DOFF = SQRT((XP0FF(I+1)-XP0FF(I))**2 + 


INT043 10 




+ ( \ POFF( 1+1) -YPOFF( I ) )**2 ) 


INT04320 




COFF(I) = (XPOFF( 1+1) -XPOFF( I) )/DOFF 


INT04330 




SOFF(I) = ( YPOFT( 1+1) -YPOFF( I) )/DOFF 


INT04340 


307 


CONTINUE 


INT04350 




GO TO 102 


INT04360 


C 


CALCULATION OF SURFACE VELOCITIES FOR THE FUNDAMENTAL SOLUTIONS 


INT04370 


C 




INT04380 


275 


DO 330 K=1 j 3 


INT04390 




DO 320 1=1, N 


INT04400 
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SUMT= T(K,I) 

SUMN=-R(K, I ) 

DO 310 J=1,N 

SUMT= SUMT + BOFF(I,J)*SIG(K,J) 

SUMN= SUMN + AOFF(I, J)*SIG(K,J) 

310 CONTINUE 

VTAN(K, I)= SUMT 
VNOR(K,I)= SUMN 
320 CONTINUE 
C 

DOFF1 = SQRT( (XPOFF( 2) -XPOFF( 1) )**2+( YPOFF( 2) -YPOFF( 1) )**2) 
DOFF 2 = SQRT( (XPOFFC 3) -XPOFF( 2) )**2+( YPOFF( 3) -YPOFF( 2) )**2) 
DOFFN = SQRT( (XPOFF(N+l) -XPOFF(N) )**2+( YPOFF(N+l) -YPOFF(N) ) 

+ v?*2) 

D0FFN1= SQRT( (XPOFF(N) -XPOFF(N-l) )**2+( YPOFF(N) -YPOFF(N-l) ) 

+ ** 2 ) 

VUT(K) = VTAN(K,N) + DOFFN * ( VTAN(K ,N) -VTAN(K,N- 1) )/ 

+ ( DOFFN+DOFFN 1 ) 

VUN(K) = VNOR(K,N) + DOFFN * (VNOR(K,N)-VNOR(K,N-l))/ 

+ ( DOFFN+DOFFN 1) 

VLT(K) = VTAN(K, 1) + DOFF1 * ( VTAN(K, 1) -VTAN(K,2) )/ 

+ (DOFF1+DOFF2) 

VLN(K) = VNOR(K, 1) + DOFF1 * ( VNOR(K, 1) -VNOR(K, 2) V 
+ (DOFF1+DOFF2) 

330 CONTINUE 

C OUTPUT FUNDAMENTAL SOLUTIONS 

IF (ICYCLE .EQ. 1 .OR. ICYCLE . GE. ICYTL-1 .OR. IP . GE. 0) 

+ WRITE( 6, 335) TITLE 
335 FORMAT! 1H1,///20A4//) 

C DO 360 K=1 , 3 

C WRITE (6,340) K 

C340 FORMAT (////, 1H /FUNDAMENTAL SOLUTION NUMBER ’,12////) 

C WRITE (6, 345) 

345 FORMAT! 3X, * I ' ,8X/X' ,11X/Y' ,10X/VT' ,10X/VN' ,8X/SIG' ///) 
C WRITE! 6,375) 1, XP(1) , YP(1), VLT(K) , VLN(K) 

C DO 350 1=1, N 

C WRITE (6, 3 75) I , X(I), Y(I), VTAN(K, I) , VNOR( K, I) , SIG(K, I ) 

C350 CONTINUE 

C WRITE! 6,375) N1 , XP(N1) , YP(N1) , VUT(K) , VUN(K) 

C360 CONTINUE 

C COMBINED FLOW AT ANGLE OF ATTACK 
C 

IF ( INLET .NE. 0) THEN 

YYY = !(VUT(3)+VLT(3))*TAN(ALPHAI)+!VUT(1)+VLT(1))*Q)/ 

+ ( ( VUT ( 3 ) +VLT( 3 ) ) - ( VUT( 2 ) +VLT( 2 ) ) *Q ) 

XXX = -( ( VUT( 1)+VLT( 1) )+( VUT( 2)+VLT( 2) )*TAN( ALPHAI) ) / 

+ ( ( VUT( 3 ) +VLT( 3 ) ) - ( VUT( 2 ) +VLT( 2 ) ) *Q ) 

ALPHA = ACOS( I. 0/SQRT! 1. 0+YYY**2) ) 

COSAL = COS! ALPHA) 

SINAL = SIN! ALPHA) 

W = XXX/SQRT(1. 0+YYY**2) 

ELSE 

COSAL = COS! ALPHA) 

SINAL = SIN! ALPHA) 

W= - ( ( VLT! 1 ) +VUT( 1 ) ) *COS AL+( VLT( 2 ) +VUT( 2 ) ) *S INAL) / 

+ (VLT! 3)+VUT( 3) ) 
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o o o 



390 

c 



395 



END IF 

FORCE COEFFICIENT CALCULATION 
SUM1= 0. 

SUMX= 0. 

SUMY= 0. 

DO 390 1=1, N 
SUM1= SUM1+ D( I ) 

SUMX= SUMX- VCOM( I )**2*S( I )*D( I ) 

SUMY= SUMY+ VCOM(I)**2*C(I)*D(I) 

CONTINUE 
FIND MAN. CHORD LENGTH 
XOMIN = XO(l) 

DO 395 I = 2 , N1 

IF ( XO(I) . GT. XOMIN) GOTO 395 

XOMIN = XO(I) 

CONTINUE 

CHORD = XO(Nl) - XOMIN 

CL1= SUM 1*25. 1327 4*W / CHORD 

CL2= ( SUMY*COSAL-SUMX*SINAL) /CHORD 

CD = ( SUMX*COSAL+SUMY*SINAL) /CHORD 

CALCULATING PARAMETERS FOR INLET VELOCITY AS MODULUS OF NOMORIZED 

IF (.NOT. SINGLE) THEN 

NUM 1 = SIN( ALPHA)+CLl*CHORD/( 4. 0*SP) 

ALPHID = ATAN2(NUM1,C0S( ALPHA)) 

NUM 1 = SIN( ALPHA) -CLl*CH0RD/(4. 0*SP) 

ALPHED = ATAN2(NUM1,C0S( ALPHA)) 

NUM1 = CLl*CHORD/( 2. 0*SP)*C0S( ALPHA) 

DNUM1= 1. 0-(CLl*CH0RD/(4. 0*SP))**2 
DALPKA = AT AN 2 ( NUM 1 , DNUM 1 ) 

UOUI = (TAN( ALPHID) -TAN( ALPHED) )*( 2. 0*SP/CHORD*COS( ALPHID) )/CLl 
CLI = CL1 * UOUI**2 
UIOU = 1. O/UOUI 

VEXIT = COS( ALPHA) /COS( ALPHED) 

ELSE 

= ALPHA 
= ALPHA 
= 0 . 0 
1. 0 
1. 0 
CLI 
1. 0 



ALPHID 
ALPHED 
DALPHA 
UOUI = 

UIOU = 

CLI = 

VEXIT = 

END IF 

FAC = 180. 0/PI 
ALPHID = ALPHID 
ALPHED = ALPHED 
DALPHA = DALPHA 
ALPHAD = ALPHA 
IF (ICYCLE .EQ. 



* FAC 

* FAC 

* FAC 

* FAC 
1 .OR. 



ICYCLE .GE. ICYTL-1 .OR. IP . GE. 0) THEN 



WRITE(6,370) ALPHAD , ALPHID, ALPHED, DALPHA, UIOU, VEXIT 



370 FORMAT (////, 1H , 'COMBINED FLOW AT 



+ 

+ 

+ 

+ 



AVERAGE 
1 



ANGLE OF 
ANGLE 



ATTACK = 
OF ' , 



F8. 3, 4X, 'DEGREES' , / , 1H ,17X,' INLET 
'ATTACK = ' ,F8. 3, 4X, 'DEGREES' ,/,lH , 

17X, 1 EXIT ANGLE = ’ ,F8. 3 ,4X, ’ DEGREES ',/, 1H ,17X, 

'TURNNING ANGLE = ' ,F8. 3 , 4X, ' DEGREES ',/, 1H ,17X, 
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365 



380 



374 



375 

C 



385 



400 



420 

430 



C 

C 

C 

C 



C 

C 

C 

C 

C 

C 

C 

C 



C 

C 

C 



+ 'INLET VEL = ' ,F10. 6,3X, 'EXIT VEL = ’ ,F10. 6,//) 

WRITE( 6 , 365 ) 

FORMAT ( 3X, ' I ' ,8X,'X0' ,10X,'Y0' ,10X,'X' ,11X,'Y' ,10X,'VT' , 

+ 10X, 'VN' ,11X, 'V' ,10X, 'CP' ,9X, 'CPI'///) 

END IF 

DO 380 1 = 1, N 

VTCOM( I )= VTAN( 1 , I )*COSAL+VTAN( 2 , I)*SINAL+W*VTAN( 3,1) 

VNCOM( I )= VNOR( 1 , I )*COSAL+VNOR( 2 , I )*SINAL+W*VNOR( 3,1) 

VCOM( I )= SQRT ( VTCOM( I ) **2 + VNC0M(I)**2) 

IF (VTCOM(I) . LT. 0.0) VCOM(I) = -VCOM(I) 

CP(I) = 1.0 - VCOM(I) ** 2 
CPI(I)= 1.0 - (VC0M(I)*U0UI)**2 

SIGCOM(I) = SIG( 1 , I)*COSAL+SIG( 2 , I)*SINAL+W*SIG( 3,1) 

XP(I) =0.5 * (X0(I)+X0(I+1)) 

YP( I) =0.5 * (Y0( I)+YO( 1+1) ) 

CONTINUE 

IF (ICYCLE . EQ. 1 .OR. ICYCLE . GE. ICYTL-1 .OR. IP . GE. 0) THEN 
WRITE (1,374) ( X0( I) , Y0( I) ,XP( I) , YP( I) ,CP(I) ,CPI(I) ,1=1, N) 

WRITE (6,375) ( I, XO(I), YO(I), XP(I), YP(I), VTCOM(I), 

+ VNCOM( I ) , VCOM( I ) , CP( I ) , CPI(I) ,I=1,N) 

FORMAT( 6F10. 4) 

FORMAT ( IX, 13, 9F12.4) 

WRITE (2) I ,X0( I) , Y0( I) ,XSTGR( I) ,YSTGR( I) ,DLS( I) ,X( I) ,Y( I) , VCOM( 
WRITE (6, 385) N+l ,X0(N+1) , Y0(N+1) 

FORMAT ( IX, 13, 2F12. 4) 

WRITE(6,400) CHORD, CL1 ,CLI 

F0RMAT(///3X,' CHORD = ' ,F10.-5,4X, 'CL(AVG) = ' ,F10. 5,4X, 

+ 'CL( INLET) = ' ,F10. 5) 

END IF 

FORMAT ( / 3X , 1HI , 6X , 2HS0 , 10X , 2HSC , 9X , 3H VN , 9X , 3HVNC , 9X , 3HDLS , 8X , 

+ 4HDLSC) 

FORMAT( 15 , 6E12. 4) 

RETURN 

END 



DATA SET KCBCAMEAN AT LEVEL 001 AS OF 08/24/84 
DATA SET KCBCAMEAN AT LEVEL 003 AS OF 04/05/84 
SUBROUTINE AMEAN( NS , ND , X , Y , IT) 



SMOOTH DATA USING 3-PTS WEIGHTING FORMULA 



NS 

ND 

X, Y 



IT 



STARTING PINT OF THE DATA TO BE SMOOTHED 
END PINT OF THE DATA TO BE SMOOTHED 
INDEPENDENT + DEPENDENT VARAIBLES OF THE DATA 
TO BE SMOOTHED 
CYCLES OF DATA SMOOTHING 



DIMENSION X(101),Y(101) 



NM = ND -NS 

IF(NM .LT. 2 .OR. IT . LT. 1) RETURN 



NDM1 = ND - 1 
NSP1 = NS + 1 
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DO 20 K=1 , IT 

DL1 = X(N'SPl) - X(NS) 

Y 1 = Y(XS) 

DO 10 I=NSP1 ,NDM1 
DL2 = X(I + 1) -X(I) 

Y2 = Y(I) 

YM = (DL2 - Yl + DL1 * Y( 1+1) )/(DLl + DL2) 

Y( I ) = 0. 5 * ( Y2 + YM) 

DL1 = DL2 

Yl = Y2 

10 CONTINUE 

20 CONTINUE 

C 

RETURN 

END 



c 


DATA 


SET 


KCBCBLGRID 


AT 


LEVEL 


001 


AS 


OF 


08/24/84 


c 


DATA 


SET 


KCBCBLGRID 


AT 


LEVEL 


001 


AS 


OF 


08/24/84 


c 


DATA 


SET 


KCBCBLGRID 


AT 


LEVEL 


004 


AS 


OF 


04/05/84 



SUBROUTINE 3LGRID(N,X,T,D1) 

C 

C GENERATE B. L. X-WISE GRID USING MODIFIED COSINE DISTRIBUTION 
C 

DIMENSION X(101),T(101),D1(101) 

DATA CRAD/57. 2957795/, BPI/3. 14159265/ 

C 

C 

c 

NN = 2 * N - 1 

EN = FLOAT((NN-l)/2) 

THO = 10. /CRAD 

CTOl = 1. + COS(THO) 

DTH = (BPI - THO) / EN 
FI = FLOAT(N - 2) 

DO 10 I=N , NN 
FI = 1. 0 + FI 

II = I - N + 1 

XII = THO + FI * DTH 

X(II) = (1.0 + COS(XII))/CTOl 
10 CONTINUE 

XI = X(l) 

XN = X(N) 

CH = XN -XI 

FN1 = FLOAT(N-l) 

N10 = N/10 

DO 20 1=1, N 

T( I ) = FLOAT( I-1)/FN1 

X(I) = (X(I)-X1)/CH 

20 CONTINUE 

C CALL SMFIT(N10 ,N,T,X,D1,N10) 

C IF(X(2). LT. 0. 35 * X(3)) X(2) = 0.35 * X(3) 

C CALL SMFIT( 1,N,T,X,D1,2) 

CALL AMEAN( 1,N,T,X,N10) 

C 

RETURN 

END 

C DATA SET KCBCBL2D AT LEVEL 001 AS OF 08/24/84 
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OOO OOO Ln ooooo 



c 

c 

c 

c 

c 



DATA SET KCBCBL2D AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCBL2D AT LEVEL 012 AS OF 04/06/84 

SUBROUTINE EL2D ( ITR , ISWPT , SURFID) 

PROGRAM CALCULATES VISCOUS/INVISCID INTERACTION USING HILBERT 
INTEGRAL. 



COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

COMMON /BLC1/ F( 101 , 2) ,U( 101, 2) , V( 101 , 2) ,W( 101 , 2) ,B( 101, 2) 

COMMON /BLC2/ DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON /BLC7/ C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 
COMMON/ EDDY 1/ RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) ,ALFAS( 100) , 
b FFS( 100) ,RTS( 100) , IEDY ,NXSPT 

COMMON /SMRY/ W( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 
b THT( 100) ,NPSTR( 100) 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /BONV/ ITMAX,EPSL,EPST,CONV 

COMMON /SAVE/ FS( 101) ,US( 101) ,VS( 101) ,WS( 101) ,BS( 101) 

COMMON /BLIN/ TITLE(20) ,XC(100) ,YC( 100) , ISG( 100) ,DELS(100) , 
b XCTR , XTR , ISTRP , ICYCLE , ICYTL,XCTRS( 2) ,TRFIND(2) 

COMMON /1SURF/ ISF 
COMMON/ PLOT /v^r v 2 ) , NXVP( 20 , 2 ) , ICC 
DIMENSION SURFID(4) 



GENERATE B. L. GRIDS + SET INITIAL CONDITIONS 

DO 5 I = 1 , NXT 
ALFAS(I) = 0. 0 
FFS(I) = 1. 0 
RTS(I) = 1. 0 
CONTINUE 

CALL INPUT( ITR, ISWPT, SURFID) 

CALCULATE HILBERT COEFFS. , C(I,J) 

CALL CALCIJCNXT, 0) 

LOOP OF CALCULATIONS 



C 

10 



20 



30 



NSS = NS 
NXSPT = NXT + 1 

IF ( ICYCLE . EQ. 1 ) NS = NXT + 1 
NX = NX + 1 

CEL =0.5 * (X(NX) + X(NX-l) ) /(X(NX) -X(NX-l)) 

CELH =0.5 * CEL 
IT =0 

RX = UE(NX)*X(NX)*RL 
SQRX = SQRT(RX) 

IT = IT + 1 

IF( IT . LE. ITMAX) GO TO 40 

NXM1 = NX-1 

CALL HEADER( TITLE , SURFID , ISTRP ) 

WRITE ( 6 , 170 ) (M,X(M),CF(M),DLS(M),UE(M),P2 (M) ,THT(M) , 

+ D(M) , ALFAS(M) , ITP(M) ,NPSTR(M) ,M=1 ,NXM1) 

WRITE ( 6 , 160 ) NX 
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STOP 






INT07210 


40 


CONTINUE 






INT07220 




IF (NX . GT. NTR ) CALL EDDY 






INT07230 




CALL COEFTR 






INTO 7240 




CALL SOLV3 






INTO 7 250 




IF( V( 1 , 2) . GT. 0. 0) GOTO 60 






INT07260 


C 








INT07270 


C 


EXTRAPOLATE CALCULATED D FOR TURBULENT SEPARATION OR LAMINAR 


INT07280 


C 


SEPARATION FOR LAMINAR FLOW CALCULATION 


ONLY 


INT07290 


c 








INT07300 




CALL EXTRAP ( NX, NXT,X,D) 






INTO 7 3 10 




NXM 1 = NX - 1 






INT07320 




CALL HEADER( TITLE , SURFID , ISTRP ) 






INT07330 




WRITE ( 6 , 170 ) (M,X(M) ,CF(M) ,DLS(M) ,UE(M) ,P2(M) ,THT(M) , 


INTO 7 340 




+ D(M) , ALFAS(M) , ITP(M) ,NPSTR(M) ,M=1,NXM1) 


INTO 7 350 




WRITE (6, 180) 






INTO 7 3 60 




WRITE ( 6 , 1 90 ) (M,X(M) ,D(M) , M=NX , NXT) 




INT07370 




GOTO 130 






INTO 7 380 


60 


IF(NX .GT. NTR) GO TO 70 






INT07390 




IF( ABS( DELV( 1) ) . GT. EPSL) GO TO 


30 




INT07400 




GO TO 80 






INTO 74 10 


70 


CONTINUE 






INT07420 




IF( ABS(DELV( 1)/V(1,2)) . GT. EPST) 


GO TO 


30 


INT07430 


80 


CONTINUE 






INT07440 


C 








INTO 7 450 


C 


CHECK FOR GROWTH 






INT07460 


C 








INT07470 




I F ( NP .GE. NP7) GO TO 90 






INT07480 




IF(ABS(V(NP,2) ) .LT. 0.0005 .AND. 


ABS( 1. 


0-U(NP-2 , 2) ) 


INT07490 




+ . LT. 0. 0035) GOTO 90 






INTO 7 5 00 




CALL FILLUP(l) 






INT07510 




IT = 1 






INT07520 




GO TO 30 






INT07530 


90 


CONTINUE 






INT07540 


C 








INT07550 




CALL FILLUP( 2) 






INTO 75 60 




CALL OUTPUT ( 1 ) 






INTO 75 70 




IF( ITR. EQ. 0 .OR. NX. GE. NTR) GOTO 


100 




INT07580 




IF(NX. LT. 3 . OR. ITR. NE. 3) GO TO 100 




INT07590 




CALL TRNS( I CODE) 






INT07600 




IF( ICODE. EQ. 1) GOTO 20 






INT07610 


100 


IF(NX .NE. NSS) GOTO 120 






INT07620 


c 








INT07630 


c 


STORE PROFILES AT THE STATION NS FOR INVERSE B. L. 


INTO 7 640 


c 


CALCULATION 






INTO 7 650 


c 








INT07660 




DO 110 J = 1 , NPT 






INT07670 




FS( J) = F( J, 2) 






INTO 7 680 




US(J) = U( J, 2) 






INT07690 




VS( J) = V( J ,2) 






INT07700 




WS(J) = W(J,2) 






INT07710 




BS( J) = B(J,2) 






INT07720 


110 


CONTINUE 






INT07730 


120 


IF ( NX . LT. NSS ) GOTO 10 






INTO 7 740 


C 


IF (ICYCLE .NE. 1 ) GO TO 130 






INTO 7 750 




IF ( NX .GE. NS ) GOTO 130 






INT07760 
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IF (NX .LT. NXT) GO TO 10 INT0777 

CALL KEADERC TITLE, SURF ID, I STRP ) INTO 7 78 

WRITE ( 6 , 170 ) (M,X(M) ,CF(M) ,DLS(M) ,UE(M) ,P2 (M),THT(M), INT077S 

+ D(M), ALFAS(M) , ITP(M) ,NPSTR(M) ,M=1 ,NXT) INT078G 

130 DO 140 I = 1 , NXT INT0781 

140 DB( I ) = D( I ) INT0782 

NS = NSS INTO 7 83 

NX = NS INTO 7 84 

NP = NPSTR(NX) INT0785 

DO 150 J = 1 , NPT INT0786 

F( J, 2) = FS( J) INT0787 

U(J,2) = US(J) INT0788 

V(J,2) = VS(J) INTO 7 89 

W(J,2) = WS(J) INT0790 

B( J, 2) = BS( J) INT0791 

150 CONTINUE INT0792 

155 INVRS = NS + 1 INTO 793 

C INT0794 

C CALCULATION SHIFTS TO USING PHYSICAL COORDINATES INT0795 

CALL MA.IN2( ITR , ISWPT , SURFID) INT0796 

C 1NT0797 

C PASS DELTA-STAR BACK TO MAIN PROG. INT0798 

C INT0799 1 

DO 158 I = 1 ,NXT INT0800' 

DELS(I) = DLS(I) INTO 8 01' 

158 CONTINUE INT0802I 

RETURN INT0803' 

C INT0804I 

C INT0805I 

C INT0806I 

160 FORMAT( 1H0, ' ** ITERATIONS EXCEEDED ITMAX AT NX = ’,15/ INT0807I 

+ 1H ,' ** CALCULATIONS STOP. ** ' ) INT0808I 

170 FORM AT ( 1H0 , ' ** SUMMARY OF STANDARD B. L. SOLUTIONS. **'/ INT0809I 

+ 1H0 , 4X , 2HNX , 7X , 1HX , 12X , 2HCF , 1 IX , 3HDLS , 12X , 2HUE , INT0810I 

+ 1 2X , 2HP2 , 1 IX , 3HTHT , 13X, 1HD, 1 OX , 4HALFA , 6X , 2HIT , 2X , 2HNP/ INT08 1 1! 

+ (1H ,3X,I3,F10. 5,2X,7E14. 5,214)) INT0812I 

180 FORMAT( 1H0 , 34H FLOW SEPARATES. D IS EXTRAPOLATED/ INT0813! 

+ 1H0 , 3X, 3H NX, 7X, 1HX, 13X, 1HD/) INT0814I 

190 FORMAT( 1H , 3X ' 13 ,F10. 5 , 2X,E14. 5 ) INT08151 

END INTO 8 161 

C DATA SET KCBCCALCIJ AT LEVEL 001 AS OF 08/24/84 INT0817I 

C DATA SET KCBCCALCIJ AT LEVEL 001 AS OF 08/24/84 INT0818I 

C DATA SET KCBCCALCIJ AT LEVEL 005 AS OF 04/05/84 INT0819' 

SUBROUTINE CALCIJ ( IL, LO) INT082O 

C INT082T 

C CALCULATE HILBERT INTEGRAL COEFFS INT0822' 

C INT0823 1 

COMMON /BLC7/ C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI INT0824' 

COMMON/EDDY 1/ RL , RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) INT0825 

+ ,ALFAS( 100) ,FFS( 100) ,RTS( 100) ,IEDY,NXSPT INT0826 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH INT0827 

DIMENSION E(2) INT0828 

C INT0829 

C INT0830 

C INT0831 1 

PI = 3. 14159265 INT0832 
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c 

c 

c 



c 

c 

c 

30 



C 

C 

c 

40 

C 

50 



60 



65 

C 

C 



C 

C 

C 

C 

C 

C 

C 



PI = PI*SQRT(RL) 

IL1 = IL - 1 
DO 65 I = 2, IL1 
E (1) = 0. 

L = LO + I 

DO 60 J = 2, IL 

J1 = J - 1 

K = J + LO 

DX1 = X(L) - X(K) 

DX2 = X(K) - X(K-l) 

DX3 = X(L) - X(K-l) 

IF ( J . EQ. I ) GO TO 30 
IF ( J .EQ. (1+1) ) GO TO 40 

J .NE. I OR 1+1 

E (2) = ( 1.0/DX2 ) * ALOGC ABS( DX3 / DX1 ) ) 

GO TO 50 

J .EQ. I 

R1 

E (2) = 

GO TO 50 

J .EQ. 1+1 



( X(K+1)-X(K) ) / (X(K+1) -X(K-l) ) 

( R1 * ALOG( ABS( DX3 / ( X(L)-X(K+1)) ) ) + 2. 0 ) / 



R1 = ( X(K- 1) -X(K-2) ) / ( X(K) -X(K-2) ) 

E (2) = ( R1 * ALOGC ABS( (X(L) -X(K-2) ) / DX1 ) ) 



2-0 ) / 



E( 2) ) / PI 



CONTINUE 

C (J1,I)= ( EC 1) 

EC 1) = E (2) 

CONTINUE 
E (2) = 0. 

J1 = IL 

K = K + 1 

C (J1,I)= EC 1) / PI 

CONTINUE 



RETURN 
END 

DATA SET KCBCCOEF 
DATA SET KCBCCOEF 
DATA SET KCBCCOEF 
SUBROUTINE COEF( GAMMA 1 ,GAMMA2) 



AT LEVEL 001 AS OF 08/24/84 
AT LEVEL 001 AS OF 08/24/84 
AT LEVEL 007 AS OF 04/05/84 



CALCULATE COEFFS OF B. L. FINITE -DIFFERENCE EQS. IN 
SEMI-TRANSF VAR;‘ T .BLES( AFTER SWITCHING). 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
COMMON /BLC1/ F < 101 , 2) ,U( 101 , 2) , V( 101 , 2) , W( 101 , 2) , B( 101 , 2) 
COMMON /BLC6/ Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
+ S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON /BLC7/ C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
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INT08440 
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INT08500 
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DX2INT08550 
INT08560 
INT08570 
INT08580 
INT08590 
INT08600 
DX2 INT08610 
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INT08650 
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OOO O --HOOO 



COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 



P1H = 0.5 * Pl(NX) 

DO 100 J= 2 ,NP 
FLARE = 1. 0 

FB = 0.5*(F(J,2) + F( J-1,2)) 

UB = 0. 5*(U( J , 2) + U( J-1,2) ) 

FVB = 0. 5*(F(J,2)*V(J,2)+F(J-1,2)*V(J-1,2)) 

IF( UB .LT. 0.0) FLARE = 0.0 
VB = 0. 5*( V( J, 2) + V( J-1,2)) 

USB = 0. 5*(U(J,2)**2 + U( J-1,2 )**2) 

WSB = 0. 5*(W(J,2)**2 + W( J-l ,2)**2) 

DERBV =(B(J,2)*V(J,2) - B( J-1,2)*V( J-1,2))/DETA( J-l) 

FB4 = 0. 5*(F( J, 1) + F( J-l , 1) ) 

VB4 = 0. 5*( V( J, 1) + V(J-1,1)) 

USB4 = 0. 5*(U(J,1)**2 + U( J-l, 1)**2) 

WSB4 = 0. 5*(W(J,1)**2 + W(J-1,1)**2) 

FVB 4 = 0. 5*(F( J, 1)*V( J, 1)+F( J-l , 1)*V( J-l , 1) ) 

DERBV4 = (B( J, 1)*V( J, 1) - B( J-l , 1)*V( J-l , 1) )/DETA( J-l) 

S1(J) = CELH*(FB - FB4) + P1H*F(J,2) + B( J,2)/DETA( J-l) 

S2( J) = CELH*(FB - FB4) + P1H*F(J-1,2) - B( J-l , 2)/DETA( J-l) 

S3( J) = CELH‘ ; '( VB + VB4) + P1H*V( J, 2) 

S4( J) = CELH- V ( VB + VB4) + P1H*V(J-1,2) 

S5( J) = -CEL*FLARE*U(J,2) 

S6( J) = - CE L- F L ARE *U (J-1,2) 

S7( J) = CEL*W( J,2) 

S8( J) = CEL‘ V W( J-l ,2) 



00 



CRB = -DERBV4 + CEL*WSB4 - CEL*FLARE*USB4 - P1(NX)*FVB4 
R2( J) = CRB - (DERBV - CEL*FLARE*USB + CEL*( VB+VB4)*(FB-FB4) 

CEL-’-WSB + Pl(NX)-'FVB) 

Rl( J) = F( J- 1 , 2) - F( J , 2) + DETA( J- 1)*UB 
R3( J-l)= U( J- 1 , 2) - U( J,2) + DETA( J-1)*VB 
R4( J- 1 )= W( J-l, 2) - W( J, 2) 

CONTINUE 



+ 



BOUNDARY CONDITIONS 

Rl( 1) = 0. 0 

R2( 1) = 0. 0 

R4(NP) = 0. 0 

I F ( NX . GE. INVRS) GO TO 120 
GAMMA 1 = 0. 0 
GAMMA2 = 1.0 
R3(NP) = 0. 0 
RETURN 
120 CONTINUE 

CII = C(NX,NX) * SQRT(X(NX) ) 

GAMMA 1 = 1. 0 

GAMMA2 = (1.0 - CII*ETA(NP))/CII 

R3( NP) = ( GI + C I I*( ETA( NP ) *W ( NP , 2 ) - F(NP,2)) -W(NP,2))/CII 
C 

RETURN 
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INT0931 
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no h-io n o noon nnnnnnn 



END 

DATA SET KCBCCOEFTR AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCCOEFTR AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCCOEFTR AT LEVEL 004 AS OF 02/21/84 

SUBROUTINE C0EFTR 

CALCULATE C0EFFS. OF B. L. FINITE -DIFFERENCE EQS. 

IN TRANSFORMED VARIABLES( BEFORE SWITCHING). 

COMMON /BLCO/ NX,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP 
COMMON /BLC 1/ F( 101 ,2) ,U( 101 ,2) ,V( 101 ,2) ,W( 101 ,2) ,B( 101 ,2) 
COMMON /BLC6/ Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
+ S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON /GRD / ETA( 101) ,DETA( 101) , A( 101) 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 



DO 100 J= 2 ,NP 

FB = 0. 5*(F( J, 2) + F( J-l ,2) ) 

UB = 0. 5*(U( J, 2) + U( J-l ,2) ) 

VB = 0. 5*(V( J,2) + V( J-1,2)) 

USB = 0. 5*(U(J,2)**2 + U( J-l ,2)**2) 

DERBV = (B(J,2)*V(J,2) - B( J-l , 2)*V( J-l , 2) )/DETA( J- 1) 
FVB = 0. 5*(F(J,2)*V(J,2) + F( J-l, 2)*V( J-1,2)) 

FVB4 = 0. 5*(F(J,1)*V(J,1) + F( J-l , 1)*V( J-l , 1) ) 

FB4 = 0. 5*(F(J,1) + F( J- 1 , 1) ) 

VB4 = 0. 5*( V( J, 1) + V(J-1,1)) 

USB4 = 0. 5*(U(J,1)**2 + U(J-1,1)**2) 

DERBV4 = (B( J, 1)*V( J, 1) - B( J-l , 1)*V( J-l , 1) ) /DETA( J-l) 



S1(J) 

S2(J) 

S3( J) 
S4( J) 
S5(J) 
S6( J) 



CELH*(FB-FB4) + 

CELH*( FB-FB4) + 
DETA(J-l) 

CELH*(VB + VB4) 

CELH*( VB + VB4) 
-(CEL+P2(NX) )*U( J,2) 

- ( C E L+P 2 ( NX ) ) *U ( J - 1 , 2 ) 



0. 5*P1(NX)*F(J,2) + B( J,2)/DETA( J-l) 
0. 5*P 1 ( NX )*F( J- 1 , 2 ) -B( J-l ,2)/ 



+ 

+ 



0. 5*P 1 ( NX )*V( J , 2 ) 

0. 5*P1(NX)*V(J-1,2) 



CLB = DERBV4 + P1(NX-1)*FVB4 - P2(NX-1)*USB4 + P2(NX-1) 

CRB = -CLB - CEL*USB4 - P2(NX) 

R2( J) = CRB - (DERBV + P1(NX)*FVB- ( CEL+P2( NX) )*USB + CEL* 
+ (VB + VB4)*(FB - FB4) ) 

Rl( J)= F( J-1,2) - F(J,2) + DETA( J-1)*UB 
R3( J-l)= U( J- 1 , 2) - U( J, 2) + DETA( J-1)*VB 
00 CONTINUE 

Rl( 1) = 0. 0 

R2( 1) = 0. 0 

R3(NP) = 0. 0 
RETURN 
END 

DATA SET KCBCCOMPBL AT LEVEL 001 AS OF 08/24/84 
DATA SET KCBCCOMPBL AT LEVEL 001 AS OF 08/24/84 
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c 

c 



c 



c 



c 



c 

c - 

c 

90 

110 

112 

120 

122 

130 

140 



150 

160 

170 

180 

182 

190 

200 

220 

210 

C 

C 

C 



DATA SET KCBCCOMPBL AT LEVEL 010 AS OF 08/24/84 



SUBROUTINE COMPBL( CASEID ,XP , YP , VT , S , DLSP ,DLS , DBPP , NBL , ITRI , XCTRI , 
+ RN ,NT, ISWPT) 



COMMON /BLCO/ NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 

COMMON /BLC7/ C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 

COMMON/ EDDY 1/ RL , RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) 

+ , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /BLIN/ TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR, ISTRP , I CYCLE , ICYTL , XCTRS( 2 ) , TRFIND( 2 ) 

COMMON /BLOW/ VN(100) 

COMMON /I SURF/ ISF 

COMMON / PLOT /NVP ( 2 ) ,NXVP( 20 , 2) , ICC 



DIMENSION 



DIMENSION 



XP( 100) ,DLSP( 100) , YP( 100) ,VT( 100) ,S( 100) , 

DBPP(IOO) ,DLS( 100) ,CASEID( 20) 

XIN( 100,2) ,YIN( 100,2) ,SI( 100 , 2) , VIN( 100 , 2) , 

DIN( 100 , 2) ,DELSTR( 100,2) ,DD( 100,2) ,DDD( 100,2) 
XB(101) ,D1( 101) ,D2(101) ,D3(101) ,IEND(2) 

NBL( 2) , ITRI( 2) ,XCTRI( 2) 



DIMENSION 
DIMENSION 
LOGICAL TRFIND 

CHARACTER * 4 SURF(4) ,STITLE(2) ,SURFID(4) 



DATA SURF / ' ' , ' R SU' , ' RFAC ' , ' E ' / 

DATA STITLE / ' UPPE ' , ' LOWE ' / 



FORMAT ( 1H1,5X, 'COMPUTING BOUNDARY LAYER USING HILBERT', 

+ ' INTEGRAL' ,/) 

FORMAT ( 1H0 ,6X, ' I ' ,9X, ' XP' , 13X, ' YP' , 13X, ' S' , 14X, ' VT' , 13X, 

+ 'DBP' / ) 

FORMAT ( 1H0,6X,'I',4X,'II',3X,'IK' ,7X,'XIN',12X,'YIN', 

+ 13X, ' SI ' , 12X, ' VIN' , 12X, 'DIN' / ) 

FORMAT ( 1H , 5X, 13 , 5E15. 6 ) 

FORMAT ( 1H , 3X, 3( 2X , 13) , 5E15. 6 ) 

FORMAT ( 1 HO, 5X, 'STAGNATION POINT IS FOUND BETWEEN POINT NO. ', 
+ 13,' AND POINT NO. ',13 / ) 

FORMAT ( 1H0,5X, 'INTERPOLATED STAGNATION POINT VALUES' / 

+ 1H0 , 5X, ' S = ' ,E13. 6 , 2X, ' XP = ' ,E13. 6 , 2X, ' YP = ', 

+ E13. 6,2X, 'DBP = ' ,E13. 6 , 2X , ' VT = 0.0' / ) 

FORMAT ( 1H0,5X, 'TOTAL NUMBER OF UPPER SURFACE POINTS = ',15, 

+ 2X, ' AND AT LOWER SURFACE = ',15 / ) 

FORMAT ( 1H0,5X, 'UPPER SURFACE DATA' ) 

FORMAT ( 1H0,5X, 'LOWER SURFACE DATA' ) 

FORMAT ( 1H0,5X, 'UPPER SURFACE CALCULATIONS' / ) 

FORMAT ( 1H0,5X, 'LOWER SURFACE CALCULATIONS' / ) 

FORMAT ( 1H0,5X, 'RESULTS OF POINT REDISTRIBUTION' / ) 

FORMAT ( 1H0,5X, 'TABLE OF DELTA-STARS' / ( 1H ,5X,8E15.6 ) ) 
FORMAT ( 1H0,5X, 'TABLE OF BLOWING-VEL. ' / ( 1H ,5X,8E15.6) ) 
FORMAT ( 1HO,5X,'NO CHANGE OF SIGN ON VT. CANNOT FIND STAG. PT. ' 
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c 






INT10560 


c 


READ ONE STRIP INPUT DATA FROM UNIT NO. 1. THE ORDER IS 


INT105 70 


c 


FROM THE LOWER SURFACE T. E. TO THE UPPER SURFACE T. E. 


INT10580 


c 






INT10590 




DO 230 I = 1,20 




INT10600 




TITLE(I) = CASEID(I) 




INT10610 


230 


CONTINUE 




INT10620 




RL = RN 




INT10630 




DO 300 I = 1 , NT 




INT10640 




DBP(I) = DBPP(I) 




INT10650 


300 


CONTINUE 




INT10660 


C 






INT10670 


C 


PRINT THE TN D U7 DATA. 




INT10680 


C 






INT10690 




P2(n = 1. 0 




INT10700 


C 


WRITE ( 6,90 ) 




INT107 10 


C 


WRITE ( 6,110 ) 




INT10720 


C 


DO 500 I = 1 , NT 




INT10730 


C 


WRITE ( 6,120 ) I,XP(I),YP(I),S(I),VT(I),DBP(I) 


INTI 0740 


C500 


CONTINUE 




INT10750 


C 






INT10760 


C 


FIND STAGNATION POINT 




INT107 7 0 


C 






INT10780 




DO 600 I = 2, NT 




INT10790 




VPROD = VT( I) * VT(I-l) 




INT10800 




IF ( VPROD . GT. 0. ) GO TO 600 




INT10810 




IS = I 




INT10820 




IS1 = IS - 1 




INT10830 




GO TO 700 




INT10840 


600 


CONTINUE 




INTI 0850 




WRITE ( 6,210 ) 




INT10860 




STOP 1 




INT10870 


C 






INT10880 


c 






INT10890 


c 


INTERPOLATE S AT VT = 0. 




INT10900 


c 






INT109 10 


c 


WRITE ( 6,130 ) IS1,IS 




INT10920 


700 


DS = S(IS) - S(IS1) 




INT10930 




DV = VT(IS) - VT(ISl) 




INT10940 




SS = SC IS) - VT(IS) * ( DS/DV ) 




INT10950 




DBB = DBP(IS) - DBP(ISl) 




INT10960 




DX = XP(IS) - XP(ISl) 




INT10970 




DY = YP(IS) - YP(ISl) 




INT10980 




DS1 = SS - S(IS) 




INT10990 




DBS = DBP(IS) + DS1 * ( DBB/DS ) 




INTI 1000 




XPS = XP(IS) + DS1 * ( DX /DS ) 




INTI 10 10 




YPS = YP(IS) + DS1 * ( DY /DS ) 




INTI 1020 


C 


WRITE ( 6,140 ) SS, XPS, YPS, DBS 




INTI 1030 


C 






INTI 1040 


C 


IU IS THE TOTAL UPPER SURFACE POINTS. + STAG. 


PT. 


INTI 1050 


C 


IL IS THE TOTAL LOWER SURFACE POINTS + STAG. 


PT. 


INTI 1060 


c 






INTI 10 70 




IU = NT - IS + 2 




INT11080 




IL = IS 




INTI 1090 


c 


WRITE ( 6,150 ) IU, IL 




INTI 1100 


c 






INTI 1110 



S3 



c 


GROUP THE DATA FOR EACH 


SURFACE. FIRST, UPPER. 


INT1112C 


c 






INT1113C 




DO 1200 L = 1,2 




INT1114C 




GO TO ( 800,900 ), L 




INT1115C 


c 






INT1116C 


c 


L = 1 IS UPPER SURFACE 




INT1117C 


c 


L = 2 IS LOWER SURFACE 




INT1118C 


c 






INT1119C 


800 


Ml = IS 




INT1120C 




M2 = NT 




INT1121C 




IEND(L) = IU 




INT1122C 




GO TO 1000 




INT1123C 


900 


Ml = 1 




INT1124C 




M2 = IL-1 




INT1125C 




IEND(L) = IL 




INT1126C 


C 






INT1127C 


1000 


1 = 1 




INT1128C 




XIN( I , L) = XPS 




INT1129C 




YIN( I , L) = YPS 




INT1130C 




SI ( I , L) = 0. 




INT1131C 




DIN( I , L) ~ DBS 




INT1132C 




VIN( I , L) = 0. 




INT1133C 




IF ( IP .GE. 1 ) THEN 




INT1134C 




IF ( L .EQ. 1 ) WRITE ( 


6,160 ) 


INT1135C 




IF ( L .EQ. 2 ) WRITE ( 


6,170 ) 


INT1136C 




WRITE ( 6,112 ) 




INT1137C 




WRITE ( 6,122 ) 1,1, I ,XIN( 1 ,L) , YIN( 1,L),SI(1,L),VIN(1,L), 


INT1138C 




+ DIN( 1 , L) 




INT1139C 




END IF 




INT1140C 




DO 1100 II = Ml, M2 




INT1141C 




I =1 + 1 




INT1142C 




IK = II 




INT1143C 




IF ( L .EQ. 2 ) IK = IL 


- II 


INT1144C 




XIN( I ,L) = XP(IK) 




INT1145C 




YIN( I , L) = YP(IK) 




INT1146C 




SI ( I , L) = S( IK) -SS 




INT1147C 




IF ( L .EQ. 2 ) SI( I ,L) 


= SS-S(IK) 


INT1148C 




VIN( I , L) = ABS( VT( IK) ) 




INT1149C 




DIN( I ,L) = DBP(IK) 




INT1150C 




IF( IP .GE. 1)WRITE ( 6, 


122 ) I , II , IK,XIN( I ,L) , YIN( I,L),SI(I,L), 


INT1151C 




+ VIN( I ,L) ,DIN( I ,L) 


INT1152C 


1100 


CONTINUE 




INT1153C 


C 






INT1154C 


1200 


CONTINUE 




INT1155C 


C 






INT1156C 


C 


RE-DISTRIBUTE POINTS ON 


EACH SURFACE TO A DENSER NUMBER. 


INT1157C 


C 






INT1158C 


C 


WRITE ( 6,90 ) 




INT1159C 




DO 2000 ISF = 1,2 




INT1160C 




NN = IEND(ISF) 




INT1161C 




ITR = ITRI(ISF) 




INT1162C 




NXT = NBL(ISF) 




INT1163C 




XCTR = XCTRI(ISF) 




INT1164C 




SURF (1) = STITLE(ISF) 




INT1165C 




ICC = 1 




INT1166C 




DO 1220 J = 1,4 




INT1167C 



S4 





SURFID(J) = SURF(J) 


INTI 1680 


1220 


CONTINUE 


INT11690 


C 


IF ( ISF .EQ. 1 ) WRITE ( 6,180 ) 


INT11700 


C 


IF ( ISF .EQ. 2 ) WRITE ( 6,182 ) 


INTI 17 10 




SCALE = SI(NN, ISF) 


INTI 1 720 


C 




INTI 17 30 




CALL BLGRID ( NXT,XB,D1,D2 ) 


INTI 1740 


C 




INTI 1750 




DO 1300 I = 1 ,NXT 


INT11760 


1300 


X (I) = XB(I) * SCALE 


INTI 17 70 


C 




INT11780 


C 


INTERPOLATE S,VT,X,Y,D INTO THE NEW DISTRIBUTION. 


INTI 17 90 


C 




INTI 1800 




CALL SMFIT ( 1 ,NN,SI( 1 , ISF) ,VIN( 1 , ISF) ,D1 ,2 ) 


INTI 18 10 




CALL SMFIT ( 1 ,NN, SI( 1 , ISF) , DIN( 1 , ISF) ,D1 , 2 ) 


INTI 1820 




CALL DIFF3 ( NN , SI( 1 , ISF) , VIN( 1 , ISF) ,D1 , D2 ,D3 , 0 ) 


INTI 1830 




CALL INTRP3( NN , SI( 1 , ISF) , VIN( 1 , ISF) ,D1 , D2 ,D3 ,NXT,X ,UE ) 


INT11840 




CALL DIFF3 ( NN , SI( 1 , ISF) , DIN( 1 , ISF) ,D1 ,D2 , D3 , 0 ) 


INTI 1850 




CALL INTRP3( NN, SI( 1 , ISF) ,DIN( 1 , ISF) ,D1 , D2 ,D3 ,NXT,X,DB ) 


INT11860 




CALL DIFF3 ( NN , SI( 1 , ISF) , XIN( 1 , ISF) ,D1 ,D2 ,D3 , 0 ) 


INTI 1870 




CALL INTRP3( NN ,SI ( 1 , ISF) , XIN( 1 , ISF) ,D1 ,D2 , D3 ,NXT,X,XC ) 


INTI 1880 




CALL DIFF3 ( NN , SI( 1 , ISF) , YIN( 1 , ISF) ,D1 , D2 , D3 , 0 ) 


INTI 1890 




CALL INTRP3( NN , SI( 1 , ISF) , YIN( 1 , ISF) , D1 , D2 ,D3 ,NXT,X , YC ) 


INTI 1900 




IF ( IP . GE. 1 ) THEN 


INTI 19 10 




WRITE ( 6,190 ) 


INTI 1920 




WRITE ( 6,110 ) 


INT11930 




DO 1320 I = 1 , NXT 


INTI 1940 




WRITE ( 6,120 ) I ,XC( I) , YC( I) ,X( I) , UE( I) ,DB( I ) 


INTI 1950 


1320 


CONTINUE 


INT11960 




END IF 


INTI 19 70 


C 




INTI 1980 


C 


INPUT TO THE B. L. PROGRAM X,UE,DB,XC, YC ARE NOW DEFINED. 


INT11990 


C 




INT12000 




DO 1350 I = 1 ,NXT 


INT120 10 




DBP(I) = DB( I ) 


INT12020 




D( I) = DB( I ) 


INT12030 


1350 


CONTINUE 


INT12040 


C 




INT12050 




CALL BL2D( I TR , I SWPT , SURF ID) 


INT12060 


C 




INT12070 




CALL DIFF3 ( NXT, X , DELS ,D1 , D2 ,D3 , 0 ) 


INT12080 




CALL INTRP3( NXT, X, DELS ,D1 , D2 ,D3 ,NN, SI( 1 , ISF) ,DELSTR( 1 , ISF) ) 


INT12090 




CALL DIFF3( NXT,X , D , D1 , D2 , D3 , 0) 


INT12100 




CALL INTRP3( NXT ,X ,D , D1 , D2,D3,NN,SI(1,ISF),DD(1,ISF)) 


INT12110 




CALL DIFF3(NN, SI( 1 , ISF) ,DD( 1 , ISF) ,DDD( 1 , ISF) , D2 , D3 , 0) 


INT12120 




TRFIND(ISF) = .FALSE. 


INT12130 




I F ( ITR .EQ. 3 .AND. NTR . LE. NXT) THEN 


INT12140 




XCTRS(ISF) = XCTR 


INT12150 




TRFIND(ISF) = .TRUE. 


INT12160 




END IF 


INT12170 


C 




INT12180 


2000 


CONTINUE 


INT12190 


C 




INT12200 


C 


PUT THE WO SURFACES DELTA -STARS BACK TO ONE STRIP 


INT122 10 


C 




INT12220 




DELSTR( 1,2) = 0. 5*( DELSTR( 2 , 1)+DELSTR( 2 , 2) ) 


INT12230 



85 



non non non onn 



SUM2 = 0. 0 

N1 = NX 1 

N2 = NXT 

DO 130 K = N1,N2 

SUM 2 = SUM 2 + C(K,NX)* (D(K) 

180 CONTINUE 
C 

N1 =2 

N2 = NX-1 

SUM1 = 0. 

DO 260 K = N1 ,N2 

SUM1 = SUM1 + C(K,NX)* (D(K) 

260 CONTINUE 

GI = UEO(NX) + SUM1 + SUM 2 

C 

RETURN 



INT1224C 
INT1225C 
INT1226C 
INT1227C 
INT1228C 
INT1229C 
INT1230C 
INT1231C 
INT1232C 
INT1233C 
INT1234C 
INT12350 
INT12360 
INT12370 
INT12380 
INT12390 
INT12400 
INT12410 
INT12420 
INT12430 
INT12440 
INT12450 
INT12460 
INT12470 
INT12480 
INT12490 
INT12500 
INTI 25 10 
INT12520 
INT12530 
INT12540 
INT12550 
INT12560 
I NT 125 70 
INT12580 
INT12590 
INT12600 
INT12610 
INT12620 
INT12630 
INT12640 
INT12650 
INT12660 
INT12670 
INT12680 
INT12690 
INT12700 
INT127 10 
INT12720 
INT12730 

DBP(K) ) INT12740 

INT12750 

C(NX,NX) * DBP(NX) INT12760 

I NT 127 70 
INT12780 



DELSTR( 1,1) = DELSTR( 1,2) 

DD( 1 , 2) =0.5 * ( DD( 2,1) + DD(2,2)) 

DD( 1 , 1 ) = DD( 1,2) 

DDD( 1,2) = 0.5 * (DDD(2,1) + DDD(2,2)) /SQRT(RL) 

DDD( 1,1) = DDD( 1,2) 

IL = IEND( 2) 

I = IL 

DO 2100 II = 2 , IL 
I = 1-1 

DLS(I) = DELSTR( 11,2) 

VN( I) = DDD(II,2)/SQRT(RL) 

DBPP(I) = DD( 11,2) 

2100 CONTINUE 
C 

I = IL-1 

IU = IEND(l) 

DO 2200 11=2, IU 
I =1 + 1 

DLS(I) = DELSTR( 11,1) 

VN( I) = DDD( II , 1)/SQRT(RL) 

DBPP(I) = DD( 11,1) 

2200 CONTINUE 

WRITE ( 6,200 ) ( DLS(I), 1=1, NT ) 

WRITE ( 6,220) (VN(I) , 1=1, NT) 

RETURN 
END 

DATA SET KCBCCOMPGI AT LEVEL 001 AS OF 08/24/84 

DATA SET KCECCOMPGI AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCCOMPGI AT LEVEL 003 AS OF 08/24/84 

SUBROUTINE COMPGI 

CALCULATE GI 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
COMMON /BLC7/ C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 



DBP(K) ) 



S6 



ooo h-*ooo t-*o o o oooooooooooo 



END 



DATA SET KCBCDIFF3 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCDIFF3 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCDIFF3 AT LEVEL 002 AS OF 04/05/84 

SUBROUTINE DIFF3 (N,X,F,FP,FPP,FPPP, IEND) 

DETERMINES THE DERIVATIVE OF THE INPUT FUNCTION AT THE INPUT PTS. 
DIMENSION X(101) ,F( 101) ,FP( 101) ,FPP( 101) ,FPPP(101) 



FIRST DERIVATIVES USING WEIGHTED ANGLES 
N1=N- 1 

DX=X( 2 ) -X( 1 ) 

DF=F( 2) -F( 1 ) 

ANG2= ATAN2(DF,DX) 

DL2=DX 

DO 10 1=2, N1 
ANG1=ANG2 
DL1=DL2 
11 = 1+1 

DX=X( II) -X( I ) 

DF=F( Il)-F(I) 

ANG2= ATAN2(DF,DX) 

DL2=DX 

ANG=( DL2-ANG 1+DL1*ANG2 ) / ( DL1+DL2 ) 

FP( I )= TAN(ANG) 

IF (I.NE.2) GO TO 10 
ANGI = ANG 
ANG1I = ANGI 
DLI = DL1 

CONTINUE 
ANGF = ANG 
ANG2F = ANG 2 
DLF = DL2 
IEND1 = IEND + 1 
GO TO (11,12,13), IEND1 

IEND = 0 , EXTRAPOLATE FOR END VALUES 

FP( 1 ) = 2. *(F( 2) -F( 1 ) ) /DLI - FP(2) 

FP(N) = 2. *(F(N) -F(N1) )/DLF - FP(N1) 

GO TO 20 

IEND = 1, DERIVATIVES ARE CONTINUOUS ACROSS ENDS 

12 ANG = (DLI*ANG2F + DLF*ANG1I) / (DLI + DLF) 

FP( 1 ) = TAN (ANG) 

FP(N) = FP( 1 ) 



INT12790 
INT12800 
INT128 10 
INT12820 
INTI 28 30 
INT12840 
INT12850 
INT12860 
INT12870 
INT12880 
INT12890 
INT12900 
INT129 10 
INT12920 
INT12930 
INT12940 
INT12950 
INT12960 
INT12970 
INT12980 
INT12990 
INT13000 
INT13010 
INT13020 
INT13030 
INT13040 
INT13050 
INT13060 
INT13070 
INT13080 
INTI 3090 
INT13100 
INTI 3 110 
INT13120 
INTI 3 130 
I NT 13 140 
INT13150 
INT13160 
INT13170 
INT13180 
INT13190 
INT13200 
INT13210 
INT13220 
INT13230 
INT13240 
INT13250 
INT13260 
INT13270 
INT13280 
INT13290 
INT13300 
INT13310 
INT13320 
INT13330 



S7 



INTI 3 34 
INT1335 
INT1336 
INT1337 
INT1338 
INT1339 
INT1340I 
INT134H 
INT13421 
INT1343I 
I NT 13441 
INT1345( 
INT1346( 
INT1347( 
INT1348( 
INT1349( 
INT1350( 
INT1351C 
INT1352C 
INT1353C 
INT1354C 
INT1355C 
INT1356C 
INT135 70 
INT13580 
INT13590 
INT136O0 
INTI 36 10 
INT13620 
INT13630 
INT13640 
INT13650 

C INT13660 

C CALCULATE EDDY VISCOSITY USING C . S. TWO -LAYERS EDDY INT13670 

C VISCOSITY FORMULA INT13680 



COMMON /BLCO/ NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP INT13690 

COMMON /BLC1/ F( 101 ,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101, 2) INT13700 

COMMON/BLC7/ C( 100,100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI INT13710 

COMMON/EDDY1/ RL , RX , S QRX , RXNTR , GMTR , GMTRS (100), INT13720 

+ ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT INT13730 

COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) INT13740 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH INT13750 

DIMENSION FINT(lOl) INT13760 

C INT13770 

C INT13780 

C INT13790 

JO=l INTI 3 8 00 

UDEL=0. 995*U(NP,2) INT13810 

DO 10 J=1 ,NP INT13820 

IF(U(J,2). LT. UDEL) JJ=J INT13830 

10 IF(U(J,2).LT. 0. 0) JO=J INT13840 

EDEL=ETA( JJ)+(ETA( JJ+1) -ETA( JJ) )/(U( JJ+1 ,2) -U( JJ,2) ) INT13850 

+ *(UDEL-U( JJ,2) ) INT13860 

DO 15 J=1,NP INT13870 

E T ADE L=ET A (J)/EDEL INT13880 



C 

C 

C 

13 



C 

C 

C 

20 



GO TO 20 

I END = 2, IF FIRST DERIVATIVE . LT. 0.0 RESET TO ZERO 
CONTINUE 

FP( 1) = 2. *(F( 2) -F( 1 ))/DLI - FP(2) 

IF (FP (1) .LT. 0.0) FP (1) =0.0 
FP(N) = 2. *(F(N) -F(N1) )/DLF - FP(N1) 

SECOND + THIRD DERIVATIVES USING CUBIC FITS 



30 



DO 30 1=2 ,N1 

11 =1-1 

12 = 1+1 
DX1 = X (II) 

DX2 = X (12) 

DF1 =2.0 * 

DF2 = 2. 0 * 

FPPP( I )= 3.0 * 

FPP (I)= DF1 - 
CONTINUE 

FPPP( 1)= FPPP (2) 
FPPP(N)= FPPP (Nl) 
FPP (1)= FPP (2 ) + 
FPP (N)= FPP (Nl) + 



- X (I) 

- X (I) 

((F (II) - F (I)) / DX1 - 
((F (12). - F (I)) / DX2 - 
(DF1 - DF2) / (DXl - DX2) 
DXl * FPPP (I) / 3.0 



FP 

FP 



(X 

(X 



( 1 ) 

(N) 



X (2 )) 
X (Nl)) 



FPPP 

FPPP 



(I)) 

(I)) 



DXl 

DX2 



(2 ) 
(Nl) 



RETURN 

END 

C DATA SET KCBCEDDY 

C DATA SET KCBCEDDY 

C DATA SET KCBCEDDY 

SUBROUTINE EDDY 



AT LEVEL 001 AS OF 08/24/84 
AT LEVEL 001 AS OF 08/24/84 
AT LEVEL 003 AS OF 04/05/84 



8S 





IF(ETADEL. GT. 1.0) ETADEL=1. 0 


INT13890 


15 


FINT( J)=l. 0/C I- 0+5. 5*ETADEL**6) 


INT13900 




CALL A ME AN ( 1 , JJ , ETA , FI NT , 2 ) 


INT13910 




RU = RL 


INT13920 




IF (IT .GT. 1) GO TO 20 


INT13930 




ALFAS(NX) = ALFAS(NX-l) 


INT13940 




FFS(NX) = FFS(NX-l) 


INT13950 




RTS ( NX) = RTS(NX-l) 


INT13960 


C 




INT13970 




GMTR = GMTRS(NX) 


INT13980 




IF (NX . LE. NS) RU = RL * UE(NX) 


INT13990 




RL2 = SQRT(RU * X(NX)) 


INT14000 




RL4 = SQRT(RL2) 


INT14010 




RL216 = 0. 16 * RL2 


INT14020 


20 


VMAX =0.5 * ( V( 1 ,2) + V(l,l)) 


INT14030 




DO 30 J=2 ,NP 


INT14040 




VM =0.5 * ( V( J,2)+V( J, 1) ) 


INT14050 




IF(VM .GT. VMAX) VMAX = VM 


INT14060 


30 


CONTINUE 


INT14070 




IF ( IEDY . EQ. 0) C-0 TO 35 


INT14080 




IF ( IT . LE. 1 . OR. GMTR . LT. 0. 85 . OR. NX . LE. NTR+3 ) 


INT14090 




+ GO TO 35 


INT14100 


C 




INT14110 


C 


MODIFY ALFA USING SIMPSON'S ARGUMENTS 


INT14120 


n 


CALL SMPSON 


INT14130 


35 


ALFA = ALFA S( NX) 


INT14140 




EDVO = ALFA * RL2 * GMTR * (U(NP,2)*ETA(NP) - F(NP,2)) 


INT14150 




DO 40 J=2 jNP 


INT14160 




JJ = J 


INT14170 




YBA = RL4*SQRT(VMAX)/26. 0*ETA(J) 


INT14180 




EL = 1. 0 


INT14190 




IF( YBA .LT. 10.0) EL = 1. 0 - EXP( -YBA) 


INT14200 




ED VI = RL216 * GMTR * ( EL*ETA( J) )**2 * ABS(V(J,2)) 


INT14210 




IF( EDVI .GT. EDVO) GO TO 70 


INT14220 




B( J,2) = 1.0 + EDVI*FINT( J) 


INT14230 




IF( B( J , 2) .LT. B( J-1,2)) B(J,2) = B(J-1,2) 


INT14240 


C 


B( J,2) = 1. 0 + EDVI 


INT14250 


40 


CONTINUE 


INT14260 




JM = 2 


INT14270 




BM = B( 2 , 2) 


INT14280 




DO 50 J=2 ,NP 


INT14290 




IF(BM. GT. B( J,2) ) GOTO 50 


INT14300 




JM = J 


INT14310 




BM = B(J,2) 


INT14320 


50 


CONTINUE 


INT14330 




GOTO 90 


INT14340 


70 


DO 80 J=JJ,NP 


INT14350 


80 


B(J,2) = 1.0 + EDVO*FINT(J) 


INT14360 


C 


80 B( J , 2) = 1.0 +EDVO 


INT14370 


C 




INT14380 


90 


CONTINUE 


INT14390 




B( 1,2) = 1. 0 


INT14400 


c" 




INT14410 




JJ = 1 


INT14420 




DO 100 J=2 ,NP 


INT14430 
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IF(U( J ,2) . LT. 0.0) JJ = J I NT 14440 

100 CONTINUE INT14450I 

IF(JJ.EQ. 1) GO TO 110 INT14460 

C INT1447C 

C IN THE SEPARATED REGION, EDDY VISCOSITY IS SET EQUAL TO INT14480 

C THAT IN THE PREVOUS STATION TO AVOID NUMERICAL TROUBLE INT14490 

JJP3 = JJ + 3 INT14500 

JJP3 = MIN( JJP3 , NP) INT145 10 

CALL AMEAN( 1 , JJP3 ,ETA,B( 1 , 2) , 2) INT14520 

110 CALL AMEAN( 1 ,NP,ETA,B( 1,2) , 1) INT14530 

C INT14540 

RETURN INT14550 

END INT14560 

C DATA SET KCBCEXTRAP AT LEVEL 001 AS OF 08/24/84 INT14570 

C DATA SET KCBCEXTRAP AT LEVEL 001 AS OF 08/24/84 INT14580 

C DATA SET KCBCEXTRAP AT LEVEL 008 AS OF 02/13/84 INT14590 

SUBROUTINE EXTRAP ( NX, NXTE,X,Y) INT14600 

C INT14610 

C EXTRAPOLATE DATA USING PARABOLIC FORMULA INT14620 

DIMENSION X( 101) , Y( 101) INT14630 

C INT14640 

Y 1 = Y(NX-2) INT14650 

Y2 = Y(NX-l) INT14660 

XI = X(NX-2) INT14670 

X2 = X(NX-l) I NT 14680 

X3 = X(NXTE) INT14690 

XI = XI -X3 INT14700 

X2 = X2 -X3 INT147 10 

BB = ( Y1-Y2) /(Xl**2 - X2**2) INT14720 

AA = Y1 - BB * Xl**2 INT14730 

DO 10 I=NX,NXTE INT14740 

XI = X(I) -X3 INT14750 

Y( I ) = AA + BB * XI** 2 INT14760 

10 CONTINUE INT14770 

RETURN INT14780 

END INT14790 

C DATA SET KCBCFILLUP AT LEVEL 001 AS OF 08/24/84 INT14800 

C DATA SET KCBCFILLUP AT LEVEL 001 AS OF 08/24/84 INT14810 

C DATA SET KCBCFILLUP AT LEVEL 007 AS OF 04/05/84 INT14820 

SUBROUTINE FILLUP( INDEX) INT14830 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS,NS , IP INT14840 

COMMON /BLC1/ F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) INT14850 

COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) INT14860 

C INT14870 

C INT14880 

C INT14890 

IF(NP. GE.NPT ) RETURN INT149O0 

IF( INDEX. EQ. 2) GOTO 10 INT14910 

C INT14920 

C DEFINE PROFILES FOR B. L. GROWTH INT14930 

NP1 = NP + 1 INT14940 

NP = NP + 2 INT14950 

NP = MIN(NP, NPT) INT14960 

NM = NP INT14970 

GOTO 20 INT14980 

C INT14990 



90 



C FILL UP PROFILES BEFORE MOVING TO THE NEXT STATION 
C 

10 NP1 = NP + 1 
NM = NPT 

20 DO 30 J=NP1 ,NM 

F(J,2) = F( J- 1 , 2 ) + DETA(J-1)*U(J-1,2) 

U(J,2) = U( J-1,2) 

V(J,2) = 0. 0 
B(J,2) = B( J-1,2) 

W(J,2) = V( J-1,2) 

30 CONTINUE 

C 

RETURN 

END 

C DATA SET KCBCHEADER AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCHEADER AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCHEADER AT LEVEL 001 AS OF 04/05/84 

SUBROUTINE HEADER ( TITLE , SURFID, ISTRP ) 

COMMON / I SURF / ISF 
C 

DIMENSION TITLE( 20) , SURFID(4) 

C 

10 FORMAT ( 1H1 , 20X, 20A4 ) 

20 FORMAT ( 1H0 , 15X, ' BOUNDARY LAYER CALCULATION FOR 

+ 'UPPER SURFACE ' , 10X , ' ICYCLE=' , 15 / 16X,71(1H-) / ) 

30 FORMAT ( 1H0 , 15X, ' BOUNDARY LAYER CALCULATION FOR 

+ ’LOWER SURFACE ' , 10X, ' ICYCLE=' , 15 / 16X,71(1H-) / ) 

C 

C 

C 

WRITE ( 6,10 ) TITLE 
IF( ISF .EQ. 1) WRITE ( 6,20 ) ISTRP 
IF( ISF .EQ. 2) WRITE ( 6,30 ) ISTRP 
C 

RETURN 

END 



c 


DATA 


SET 


KCBC INPUT 


AT 


LEVEL 


001 


AS 


OF 


08/24/84 


c 


DATA 


SET 


KCBC INPUT 


AT 


LEVEL 


001 


AS 


OF 


08/24/84 


c 


DATA 


SET 


KCBC INPUT 


AT 


LEVEL 


009 


AS 


OF 


08/24/84 



SUBROUTINE INPUT( ITR , ISWPT, SURFID) 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS,NS, IP 

COMMON /BLC1/ F( 101 ,2) ,U( 101 ,2) , V( 101 , 2) ,W( 101 ,2) ,B( 101 , 2) 

COMMON /BLC2/ DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON /BLC7/ C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
COMMON /BONV/ ITMAX,EPSL,EPST,CONV 
C0MM0N/EDDY1/ RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) , 

+ ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /BLIN/ TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR , ISTRP , I CYCLE , ICYTL, XCTRS( 2 ) ,TRFIND( 2 ) 

COMMON/TRN/ PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON /I SURF/ ISF 
DIMENSION D 1 ( 100) ,D2( 100) ,D3( 100) 

DIMENSION SURFID(4) ,XCS( 100) 

LOGICAL TRFIND 
C 
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ITMAX 


— 


15 


EPSL 




0. 0001 


EPST 


= 


0. 01 


NPT 


= 


101 


ETAE 


= 


8. 0 


VGP 


= 


1. 14 


DETA1 


= 


0. 01 


NSS 


= 


NXT / 



C SEARCH FOR PRESSURE PEAK AS SWITCH POINT 
UMAX = UE( 1 ) 

DO 50 I = 2 , NXT 

IF ( UE ( I ) .LE. UMAX) GO TO 55 

UMAX = UE( I ) 

50 CONTINUE 

GO TO 60 

55 NS = I - 4 

60 IF (NS .GT. NSS) NS = NSS 

IF (NS .LT. 3) NS = 3 

C 

C CALCULATE THE PRESSURE PARAMETERS PI + 

C USING TRANSFORMED COORDINATES 
C 

CALL DIFF3 (NXT, X, UE, Dl, D2, D3, 0 
DO 65 I = 2, NXT 
P2( I ) = X(I) * Dl( I ) /UE( I) 

P1(I) = 0. 5 * (1. 0 + P2( I ) ) 

65 CONTINUE 

Pl(l) = 0. 5 * (1. 0 + P2( 1 ) ) 

XCMIN = XC(1) 

MIN = 1 

DO 70 1=1, NXT 
IF(XCMIN. LT. XC( I ) ) GOTO 70 
XCMIN = XC(I) 

MIN = I 

70 CONTINUE 

DO 80 I = 1 , NXT 
IF (I . GE. MIN) THEN 
XCS(I) = XC(I) 

ISG(I) = 1 
ELSE 

XCS(I) = -XCS(I) 

ISG(I) = -1 
END IF 

80 CONTINUE 

INVRS = NS + 1 
C 

C SEARCH FOR TRANSITION LOCATION 
ITRP1 = ITR + 1 
GOTO (150, 95, 120, 150), ITRP1 
C 

C TRANSITON LOCATION IS INPUT ( = XCTR) 

95 DO 100 1=1, NXT 

IF(XCTR. LT.XCS(I)) GOTO 105 
100 CONTINUE 
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NTR = NXT + 1 


INT16120 




GOTO 200 


INTI 6 130 


105 


NTR =1-1 


INT16140 




IF (NTR. LT. 3) THEN 


INT16150 




NTR = 3 


INT16160 




XCTR = XC(NTR) 


INT16170 




END IF 
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GOTO 200 


INT16190 


C 
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c 


TRANSITION LOCATION IS SET AT THE PRESSURE PEAK 


INT16210 


c 
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120 


UN = UE( 1) 


INT16230 




IN = 1 


INTI 6 240 




DO 75 I = 1 ,NXT 
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IF(UM. GT.UE(I)) GOTO 75 
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IN = I 


INT16270 




UN = UE(IN) 
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75 


CONTINUE 
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IF( IN. LT. 4) IN = 4 
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NTR = IN 


INT16310 




XCTR = XC(NTR) 
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GOTO 200 
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C 
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TRANSTION LOCATION WILL BE CALCULATED BASED ON NICHEL CRITERION 


INT16350 


C 




INT16360 


150 


NTR = NXT + 1 


INT16370 


200 


DO 210 1=1, NXT 
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210 


GNTRS( I)= 0.0 


INT16390 


C 




INT16400 


C ' 


TRANSITION LOCATION PROVISIONALLY FROM PREVIOUS CYCLE 


INT16410 


C 




INT16420 




IF (TRFIND(ISI) ) THEN 


INT16430 




DO 211 I = 1 , NXT 
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XCS(I) = XO'T) 


INT16450 




IF (I .Li. NIN) XCS(I) = -XCS(I) 


INT16460 


211 


CONTINUE 


INT16470 




DO 215 1=1, NXT 


INTI 6480 




IF(XCTRS( ISF) .LE.XCS(I)) GOTO 217 


INT16490 


215 


CONTINUE 


INT16500 


217 


NTR =1-1 


INT16510 




XCTR = XCTRS(ISF) 
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IF (NTR .LT. 3) THEN 
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NTR = 3 


INT16540 




XCTR = XC(NTR) 
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END IF 


INT16560 




END IF 


INT16570 


C 




INT16580 


C 


CALCULATE GAMTR DISTRIBUTION 


INTI 65 90 


c 




INT16600 




IF (NTR. GT.NXT-1) GOTO 250 


INT16610 




FAC = (XCTR-XC(NTR) )/ (XC(NTR+1) -XC(NTR) ) 


INT16620 




XTR = X(NTR) + FAC*(X(NTR+1) -X(NTR) ) 


INT16630 




UETR = UE(NTR) + FAC*(UE(NTR+1) -UE(NTR) ) 


INT16640 




RXNTR = XTR * UETR * RL 


INT16650 




GGFT = 1. 0/PGAMTR*RL**2/RXNTR**l. 34 
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DO 220 I=NTR+1 ,NXT 
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ALFAS(I) = 0.0168 
220 GMTRS( I)= 1.0 

ALFAS(NTR) = 0.0168 
UEINTG = 0. 0 
U1 = 0.5/UETR 
XI = XTR 

DO 230 I=NTR+1 ,NXT 
U2 = 0. 5/UE( I ) 

X2 = X(I) 

UEINTG = UEINTG+(U1+U2)*(X2-X1) 

U1 = U2 
XI = X2 

GG = GGFT*UE I NTG* ( X ( I ) -XTR)*UE( I )**3 
IF(GG . GT. 10. 0) GO TO 250 
GMTRS(I) = l.O-EXP(-GG) 

230 CONTINUE 
250 CONTINUE 
C 

C GENERATE B. L. ETA GRIDS + INTIAL VELOCITY PROFILES 
C 

CALL INTL(ETAE,DETA1 ,VGP) 

DO 260 1=1, NXT 
UEO(I) = UE( I) 

260 CONTINUE 
C PRINT OUT INPUT DATA 
C 

IF (ICYCLE .GE. ICYTL-1 .OR. IP . GE. 0) THEN 
CALL HEADER( TITLE , SURFID , ISTRP ) 

WRITE( 6 , 1002) NXT, ITR, IP, NS ,NTR, ISWPT 
WRITE ( 6 , 1 00 3 ) VGP, DETA1 , RL , XCTR , OMEGA , PGAMTR 
ELSE 

IF (ISF.EQ. 1) WRITE( 6 , 1008) ICYCLE 
IF (ISF.EQ. 2) WRITE( 6 , 1009) ICYCLE 
END IF 

IF (NTR. LT. NXT) THEN 

IF (ITR.EQ. 1) WRITE (6,1005) XCTR, XTR, NTR 
IF (ITR. EQ. 2) WRITE (6,1006) XCTR, XTR, NTR 
IF (TRFIND(ISF)) WRITE(6, 1007) XCTR, XTR, NTR 
END IF 
RETURN 



C 

C 

c 

2 

3 



FORMAT( 20A4) 
FORMAT(10L c n 



4 

1001 

1002 


FOR.MAT(6F10 
FORMAT( 1H1 , 
FORMAT ( 1H0 , 


.0) 

20X, 

10H 


20A4) 

NXT 


= , 15 , 7X, 10H ITR 


= ,15, 


7X/ 




+ 


1H 


, 10H 


IP = , 15 , 7X, 10H 


NS 


, 15 , 7X/ 




+ 


1H 


, 10H 


NTR = , 15 , 7X, 10H 


I SWPT= 


,15) 


1003 


FORMAT ( 1H0 , 


10H 


VGP 


= ,E12. 4, 10H DETA1= ,E12. 4/ 




+ 


1H 


, 10H 


RL = ,E12. 4, 10H 


XCTR = 


,E12. 4/ 




+ 


1H 


, 10H 


OMEGA = ,E12. 4, 10H 


PGAMTR= 


,E12. 4) 



1004 FORMAT( 1H0 , 3X , 2H I , 6X , 2HXC , 1 IX , 2HYC , 1 IX , 2H X , 1 IX , 2HUE , 1 IX , 2HP 1 , 

+ 1 IX, 2HP2 , / ( 1H , 3X, 13 , 6E13. 5 ) ) 

1005 FORMAT(/3X, 'BEGIN OF TRANSITION IS BEING INPUT AT X/C =',F8.4,4X 



INT16680 
INT16690 
INT16700 
INT16710 
INT16720 
INT16730 
INT16740 
INT16750 
INT16760 
INT167 7 0 
INT16780 
INT16790 
INT16800 
INT168 10 
INT16820 
INT16830 
INT16840 
INT16850 
INT16860 
INT16870 
INT16880 
INT16890 
INT16900 
INT169 10 
INT16920 
INT16930 
INT16940 
INT16950 
INT16960 
INT16970 
INT16980 
INT16990 
INT17000 
INT17010 
INT17020 
INT17030 
INT17040 
INT17050 
INT17060 
INT17070 
INT17080 
INT17090 
INT17100 
INT17110 
INT17120 
. INT17130 
INT17140 
INT17150 
INT17160 
INTI 7 170 
INT17180 
INT17190 
INT17200 
INT17210 
INT17220 
, INT17230 



94 



+ 'S/C =' ,F8. 4,4X, 'NTR =' ,13/) 

1006 FORMAT(/3X, 'BEGIN OF TRANSITION IS SET AT PRESSURE PEAK, X/C =' , 

+ F 8 . 4 , 4X , ' S / C =' ,F8. 4,4X, 'NTR =',I3/) 

1007 FORMATC/3X, 'BEGIN OF TRANSITION IS PROVISIONALLY TAKEN FROM ', 

+ 'PREVIOUS CYCLE: X/C=' ,F8. 4,4X, ' S/C =' ,F8. 4,4X, 'NTR =',I3/) 

1008 FORMATC//3X, 'UPPER SURFACE CALCULATIONS OF CYCLE', 13) 

1009 FORMAT(//3X, 'LOWER SURFACE CALCULATIONS OF CYCLE', 13) 

END 

C DATA SET KCBCINTL AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCINTL AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCINTL AT LEVEL 009 AS OF 02/22/84 

SUBROUTINE INTL(ETAE,DETA1,VGP) 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

COMMON /BLC1/ F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON /BLC2/ DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON /BLC6/ Sl(101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101), 

+ S7(101),S8(101),R1(101),R2(101),R3(101),R4(101) 

COMMON /BONV/ ITMAX,EPSL,EPST,CONV 
COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
C 

C 

c 

C GENERATE THE GRID 

C 

DETA(l) = DETA1 

IF( VGP. LT. 1. 0) VGP = 1. 0 

IF C ( VGP-1. 0) . LE. 0.001) GO TO 10 

NP = ALOG( (ETAE/DETA( 1) )*( VGP-1. 0)+l. 0)/AL0G( VGP) + 1.001 
GO TO 20 

10 NP = ETAE/DETA( 1) + 1.001 

20 IFCNP . LE. NPT) GO TO 30 

WRITE(6, 150 ) 

STOP 

30 ETA(l) = 0. 0 
DO 40 J=2 ,NPT 

ETA(J) = ETA(J-l) + DETA(J-l) 

DETA( J)= VGP‘ A 'DETA( J-l) 

A(J) = 0. 5*DETA( J-l) 

40 CONTINUE 

C 

C GENERATE INITIAL VELOCITY PROFILE 
80 DO 90 J=1 ,NP 

ETAB = ETA( J)/ETA(NP) 

ETAB2 = ETAB**2 

F( J , 2) = 0. 25*ETA(NP)*ETAB2*(3. 0 - 0. 5*ETAB2) 

U( J , 2 ) = 0. 5*ETAB*(3. 0 - ETAB2) 

V(J,2) = 1. 5*( 1. 0 - ETAB2)/ETA(NP) 

B( J , 2) = 1. 0 
W(J,2) = 1. 0 
90 CONTINUE 

NX =1 

IT =0 

120 IT = IT + 1 

IF( IT .LT. UMAX) GO TO 130 



INTI 7 240 
INT17250 
INTI 7 260 
INT17270 
INTI 7 280 
INT17290 
INT17300 
INT173 10 
INT17320 
INT17330 
INT17340 
INT17350 
INT17360 
INT17370 
INT17380 
INT17390 
INT17400 
INT17410 
INT17420 
INT17430 
I NT 17 440 
INT17450 
INT17460 
INT17470 
INT17480 
INT17490 
INT17500 
INT17510 
INT17520 
INT17530 
INTI 7540 
INT17550 
INTI 75 60 
INT17570 
INTI 7 580 
INT17590 
INT17600 
INT17610 
INT17620 
INT17630 
INT17640 
INTI 7 650 
INT17660 
INT17670 
INT17680 
INT17690 
INT17700 
INT17710 
INT17720 
INT17730 
INT17740 
INT17750 
INT17760 
INT17770 
INT17780 



95 



ooooooooooooo ooo 



STOP 

130 CONTINUE 

c 

DO 140 J= 2 ,NP 

FB = 0. 5*(F( J,2) + F( J- 1 , 2) ) 

UB = 0. 5*(U( J,2) + U( J-1,2)) 

VB = 0. 5*(V( J, 2) + V( J-1,2)) 

USB = 0. 5*(U(J,2)**2 + U( J- 1 , 2)**2 ) 

DERBV = (B(J,2)*V(J,2) - B( J-l , 2)*V( J-l , 2) )/DETA( J- 1) 
FVB = 0. 5*(F(J,2)*V(J,2) + F( J-1,2)*V(J-1,2)) 

C 

S1(J) = 0. 5*P1(NX)*F(J,2) + B( J,2)/DETA( J-l) 

S2( J) = 0. 5*P1(NX)*F(J-1,2) - B( J-l , 2)/DETA( J-l) 

S3( J) = 0. 5*P1(NX)*V(J,2) 

S4( J) = 0. 5*P1(NX)*V(J-1,2) 

S5( J) = -P2(NX)*U( J,2) 

S6( J) = -P2(NX)*U( J-1,2) 

CRB = -P2(NX) 

R2( J) = CRB - (DERBV + P1(NX)*FVB - P2(NX)*USB) 

C 

Rl( J) = F( J-l , 2) - F(J,2) + DETA( J-l) "UB 
R3(J-1)= U ( J - 1 , 2 ) - U( J,2) + DETA(J-1)*VB 
140 CONTINUE 

Rl( 1) = 0. 0 

R2( 1) = 0. 0 

R3(NP) = 0. 0 
CALL SOLV3 

IF( ABS( DELV( 1 ) ) . GT. EPSL ) GO TO 120 
CALL FILLUP( 2) 

CALL OUTPUT(l) 

C 

RETURN 

C 

150 FORMAT ( 1H0 , 37HNP EXCEEDED NPT -- PROGRAM TERMINATED) 

END 

DATA SET KCBCINTRP3 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCINTRP3 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCINTRP3 AT LEVEL 003 AS OF 04/05/84 

SUBROUTINE INTRP3 (N1 ,X1 ,F1 ,FP1 ,FPP1 ,FPPP1 ,N2 ,X2 ,F2) 

CUBIC INTERPOLATION 

GIVEN THE VALUES OF A FUNCTION (FI) AND ITS DERIVATIVES 
AT N1 VALUES OF THE INDEPENDENT VARIABLE (XI) 

FIND THE VALUES OF THE FUNCTION (F2) 

AT N2 VALUES OF THE INDEPENDENT VARIABLE (X2) 

X2 CAN BE IN ARBITRARY ORDER 



DIMENSION Xl(101) ,F1( 101) ,FP1( 101) ,FPP1( 101) ,FPPP1( 101) ,X2( 101) , 
+ F2( 101) 

DATA EPS /IE -04/ 
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20 

30 



10 

c 



JT = 2 

DO 10 I = 1 , N2 

XM = X2( I ) 

DO 20 J = JT , N 1 

J1 = J -1 

IF (Xl(J).GE.XM ) GO TO 30 
CONTINUE 
J = N1 
JT = J 

DXX = X2( I ) - X1(J1) 

DXX2 = DXX * DXX / 2. 

DXX3 = DXX2 * DXX / 3. 

F2( I ) = F1(J1) + DXX*FP 1 ( J 1 ) + DXX2*FPP1(J1) + DXX3*FPPP1( Jl) 



RETURN 

END 

DATA SET KCBCJOIN AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCJOIN AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCJOIN AT LEVEL 012 AS OF 02/20/84 

SUBROUTINE J0IN( INDEX) 

COMMON /BLCO/ NX,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP 

F(101,2),U(101,2),V(101,2),W(101,2) ,B( 101,2) 

C( 100, 100) ,D( 100) ,DB( 100) ,DBP'( 100) ,UE0( 100) ,GI 
RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) 

, ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 
COMMON /GTY / X(101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) 

COMMON /SAVE/ FS( 101) ,US( 101) ,VS( 101) ,BS( 101) ,WS( 101) 

COMMON /SMRY/ VW( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 
b THT( 100) ,NPSTR( 100) 



COMMON /BLC1/ 
COMMON /BLC7/ 
COMMON/ EDDY 1/ 



INDEX = 1 FOR THE FIRST SWEEP 
INDEX = 2 FOR SUBSEQUENT SWEEP 

CALL COMPGI 
Cl I = C(NX,NX) 

UES = GI / (1.0 - DLS(NX) * SQRT(RL) * CII ) 
IF( INDEX. EQ. 1) GOTO 15 
C 

C RETRIEVE PROFILES AT STATION NS FOR INVERSE B. L. 
C CALCULATION 

DO 10 J=1 ,NPT 
F(J,2) = FS( J) 

U(J,2) = US(J) 

V(J,2) = VS( J) 

W(J,2) = WS(J) 

B( J, 2) = BS( J) 

10 CONTINUE 

UES = UES/W( 1,2) 

SQS = 1. 0 

GOTO 30 
15 CONTINUE 

SQS = 1.0 / SQRT(UES) 

DO 20 J=2 ,NPT 
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ETA(J) = ETA( T )*SQS 
DETA(J-l) = ExA(J) - ETA(J-l) 

A(J) = 0. 5*DETA(J-1) 

20 CONTINUE 
C 

30 DO 60 J=1 ,NPT 

U(J,2) = U(J,2)*UES 
W( J , 2) = UES * W(J,2) 

F(J,2) = F( J , 2 ) *SQS*UES 
V( J , 2) = V( J, 2)/SQS*UES 
60 CONTINUE 

UE(NX) = W(l,2) 

RX = UE ( NX ) *X( NX ) *RL 
SQRX = SQRT(RX) 

C IF(NX. GT.NTR) CALL EDDY 

CALL FILLUP( 2) 

IF( INDEX. EQ. 2) GOTO 70 
C 

C STORE PROFILES AT STATION NS FOR THE NEXT SWEEP 
DO 65 J=1 ,NPT 
FS( J) = F( J, 2) 

US(J) = U( J, 2) 

VS( J) = V( J, 2) 

WS(J) = W(J,2) 

BS( J) = B( J,2) 

65 CONTINUE 

70 DO 80 J=1,NPT 

F( J, 1) = F( J, 2) 

U( J, 1) = U(J,2) 

V ( J , 1 ) = V ( J , 2 ) 

W( J, 1) = W ( J , 2 ) 

B( J , 1 ) = B( J, 2) 

80 CONTINUE 

RETURN 
END 

C DATA SET KCBCMAIN AT LEVEL 005 AS OF 09/18/84 

C 

C PROGRAM MAIN 

C 

C 
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SUBROUTINE C ASBLP( K2 , XP , YP , XMP , YMP , XS , YS , DLSP , VC , DBPP 
+ , RN,NBL, ITRI ,XCTRI , CASE ID) 

COMMON /BLCO/ NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 
COMMON/BLIN/ TITLE(20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR,XTR, ISTRP, ICYCLE , ICYTL,XCTRS(2) ,TRFIND( 2) 

COMMON/ EDDY 1/RL , RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) , 

+ ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY ,NXSPT 

COMMON/ BLOW/ VN(IOO) 

COMMON/TRN/ PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON/PLOT/NVP( 2 ) , NXVP( 20 , 2 ) , ICC 



DIMENSION 


CASEID( 20 ), 


XCTRI 


( 2 ), 


ITRI 


( 2 ) 


DIMENSION 


XP 


( 100), 


YP 


( 100), 


XMP 


( 100) 


DIMENSION 


YMP 


( 100), 


VC 


( 100), 


SM 


( 100) 


DIMENSION 


XS 


( 100), 


YS 


( 100), 


NBL 


( 2 ) 


DIMENSION 


VT 


( 100), 


S 


( 100), 


DLSP 


( 100) 


DIMENSION 


DLS 


( 100), 


XO 


( 100), 


YO 


( 100) 
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c 


DIMENSION 

DIMENSION 


! D1 ( 100), D2 ( 100), D3 ( 100) 

I XPS (100), YPS(IOO) , DBPP(IOO) 


INT19470 

INT19480 

INT19490 


10 


FORMAT 


( 


20A4 ) 


INT19500 


20 


FORMAT 


( 


415 ) 


INT195 10 


30 


FORMAT 


( 


315 , 3F10. 0 ) 


INT19520 


40 


FORMAT 


( 


6F10.5 ) 


INT19530 


100 


FORMAT 


( 


1H1 ,5X, 20A4 ) 


INT19540 


110 


FORMAT 


( 


1H0,6X, 'CYCLE NO. = ',15 ) 


INT19550 


120 


FORMAT 

+ 


( 


1H0,6X, 'FINISHED TOTAL NUMBER OF CYCLES = ',15, 
4X, ' JOB COMPLETED. ' ) 


INT19560 

INT19570 


130 


FORMAT 

+ 


( 


1H0 , 6X , ' THE UPDATED DISPLACEMENT SURFACE ',/ IX , 2HNX , 
9X , 2HXP , 12X , 2HYP , 1 IX , 3HDLS , 10X , 4HDBPP , 12X , 2HVN) 


INT19580 

INT19590 


140 


FORMAT 


( 


15 , 5E14. 6) 


INT19600 


150 


FORMAT 

+ 


( 


1H0 . 6X, ' READ IN CONTROL POINTS DATA' ,/lX,2HNX,9X, 
3HXMP , 11X, 3HYMP, 12X, 2HSM, 13X, 2HVC) 


INT19610 

INT19620 


160 


FORMAT 


( 


15 , 4E14. 6) 

1H0,6X, ' INTERPOLATED ORIGINAL GEOMETRY DATA’ ,/ IX, 
2HNX , 9X , 2HXP , 12X , 2HYP , 12X , 1HS , 14X , 2HVT) 


INT19630 


170 


FORMAT 

+ 


( 


INT19640 

INT19650 


180 


FORMAT 


( 15 ,4E14. 6) 


INT19660 



c 

C 

c 

GRANG(X1,X2,X3,Y1,Y2,Y3,X0)= (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 
+ +(X0-X1)*(X0-X3)/(X2-X1)/(X2-X3)*Y2+(X0-X1)*(X0-X2) 

+ / ( X3 -XI ) / ( X3 -X2 ) *Y3 

ISWPT = 1 
IEDY = 1 
KBL(l) = 91 
NBL( 2) = 71 

WRITE ( 6,100 ) CASEID 
CONTINUE 
ISTRP = ICYCLE 
NN = K2 - 1 
NT = K2 

WRITE ( 6,110 ) ICYCLE 

INTERPOLATE OUTPUT CONTROL POINTS TO ORIGINAL GEOMETRY 

DO 15 I = 1 , NT 
XPS(I) = XP(I) 

YPS(I) = YP( I ) 

15 CONTINUE 

S( 1) = 0. 0 
DO 50 I = 2 , NT 

S( I) = S(I-l) + SQRT( (XP( I) -XP( 1-1) )**2 + 

+ ( YP( I) -YP( 1-1) )**2) 

50 CONTINUE 

SM(1) = SQRT( (XMP( 1) -XP( 1) )**2 + ( YMP( 1) -YP( 1) )**2) 

DO 60 I = 2 , NN 

SM( I ) = SM(I-l) + SQRT((XMP(I)-XMP(I-1))**2+ 

+ ( YMP( I) -YMP( 1-1) )**2) 

60 CONTINUE 

C CALL AMEAN(NN-10,NN,SM,VC, 1) 

CALL AME AN( 1 , NN , SM , VC , 1 ) 

SNT = S ( NT ) 



C 

C 

5 



C 

C 

c 
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SM(NT) = SM(NN)+SQRT( (XMP(NN) -XP(NT) )**2+( YMP(NN) -YP(NT) )**2) 
SMNT = S ( NT ) / SM(NT) 

DO 65 I = 1 , NN 
65 SM( I ) = SM( I ) * SMNT 

C ALL DIFF3 ( NN , SM , VC , D 1 , D2 , D3 , 0 ) 

CALL INTRP3 ( NN , SM , VC , D 1 , D2 , D3 , NT , S , VT) 

C PRINT OUT INPUT DATA 
C 

C WRITE (6 , 150) 

C WRITE( 6 , 160) ( I ,XMP( I) , YMP( I) ,SM( I) , VC( I) , 1=1 ,NN) 

C WRITE( 6, 170) 

C WRITE( 6 , 180) ( I ,XP( I) ,YP( I) ,S(I) ,VT(I) ,1=1, NT) 

C 

XPMIN = XP(1) 

DO 44 I = 2 , NT 

IF (XP(I) . GT. XPMIN) GO TO 44 

XPMIN = XP(I) 

44 CONTINUE 

CHORD = XP(NT) - XPMIN 

DO 45 I = 1 , NT 

XP(I) = (XP( I ) -XPMIN) / CHORD 

YP( I) = YP(I) / CHORD 

45 CONTINUE 

CALL COMPBL ( CASE ID , XP , YP , VT , S , DLSP , DLS , DBPP , NBL , ITRI , XCTRI , 
+ RN , NT , I SWPT ) 

C 

C 

C SMOOTH THE CALCULATED DISPLACEMENT THINKNESS 
C 

C CALL SMFIT( 1 , NT , S , DLS , D 1 , 2 ) 

C 

C ADD DISPLACEMENT THINKNESS ON THE ORIGINAL BODY 
C 

DO 70 I = 1 , NT 
DLSP(I) = DLS(I) 

70 CONTINUE 

C CALL SMFIT ( 1 , NT,XS , YS ,D1 , 2) 

DO 80 I = 1 , NT 
XP(I) = XPS(I) 

YP( I) = YPS(I) 

80 CONTINUE 

IF (ICYCLE . GE. ICYTL-1 .OR. IP . GE. 0) THEN 

C WRITE (6,130) 

C WRITE (6,140) ( I ,XP( I) , YP( I) ,DLS( I) ,DBPP( I) , VN( I) , 1=1 ,NT) 

WRITE ( 6,120 ) ICYCLE 
END IF 
C 

RETURN 

END 

C DATA SET KCBCMAIN2 AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCMAIN2 AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCMAIN2 AT LEVEL 010 AS OF 04/06/84 

SUBROUTINE MAIN2( ITR, ISWPT, SURFID) 

COMMON /BLCO/ NX , NXT , NP , NPT , NTR , IT , INVRS ,NS , IP 

COMMON /BLC1/ F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON /BLC2/ DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
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c 

c 



10 

15 

20 



25 



30 

70 

C 



COMMON /BLC7 / C( 100,100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
COMMON /BONY/ ITMAX , EPSL , EPST , CONY 
COMMON/EDDY 1 / RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) 

+ , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /SMRY/ W( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 

+ THT( 100) ,NPSTR( 100) 

COMMON /BLIN/ TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR , ISTRP , ICYCLE , ICYTL,XCTRS( 2) ,TRFIND( 2) 

COMMON/ BLC 9/ UEB(IOO) , CFS(IOO) 

COMMON/TRN/ PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON /I SURF/ ISF 
COMMON/PLOT/NVP( 2) ,NXVP( 20 , 2) , ICC 
DIMENSION SURFID(4) ,RTSS( 11) 

LOGICAL SMOOTH , SEPART , HOMOPY 

DATA RTSS/O. 0 , 0. 1 , 0. 2 , 0. 3 , 0. 4 , 0. 5 , 0. 6 , 0. 7 ,0. 8 , 0. 9 , 1. 0/ 



GRANG(X1 ,X2 ,X3,Y1,Y2,Y3,X0)= (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 
+ +( XO -X 1 )*( XO -X3) /( X2 -X 1 ) / ( X2 -X3 ) *Y2+( XO -X 1 ) *( XO -X2 ) 

+ /(X3-X1)/ (X3-X2)*Y3 

ISWP = 0 

INDEX = 1 

IGROWT = 2 

NXSPT = NXT + 1 
CALL JOINC INDEX) 

NXSTOP = NXT-1 
IF (NS . GE. NTR ) GOTO 15 
ISWP = ISWP + 1 

NX = NX + 1 

HOMOPY = PAT 

CEL = 0. 5*(X(NX)+X(NX-1))/(X(NX)-X(NX-1)) 

PI (NX) =0.5 
P2(NX) = 0. 0 

CELH = 0. 5*CEL 

IT =0 

CALL COMPGI 
IGR0W=1 

IT = IT + 1 

RX = UE ( NX ) *X ( NX ) *RL 

SQRX = SQRT(RX) 

IF( IT . LE. ITMAX) GO TO 80 
IF( HOMOPY) GO TO 72 
IRC = 1 

RT = RTSS (IRC) 

HOMOPY = .TRUE. 

UEREF = UEO(NX-l) 

UESAVE = UEO(NX) 

UEO(NX) = RT*UESAVE+(1.0-RT)*UEREF 
DO 61 J=1 ,NP 
F(J,2) = F( J , 1 ) 

U(J,2) = U( J , 1) 

V( J , 2) = V( J, 1) 

W ( J , 2 ) = W ( J , 1 ) 
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B(J,2) = B(J,1) 

61 CONTINUE 

GO TO 30 
C 

72 NXSTOP = NX - 1 

CALL AME AN ( NS , NXSTOP , X , CF , 1 ) 

C CALL AME AN( NS, NXSTOP, X,VW,1) 

C CALL HEADER( TITLE , SURFID , ISTRP ) 

WRITE (6, 250 ) ISWP 

WRITE (6, 260 ) (M,XC(M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 

+ UEO(M) ,D(M) ,DB(M) ,GMTRS(M) , ITP(M) ,NPSTR(M) ,M=1 , NXSTOP) 

WRITE( 6 , 270 ) NX 
STOP 

80 CONTINUE 

IF(NX . GT. NTR) GOTO 100 
C 

C LAMINAR FLOW CALCULATION 
C 

CALL COEF( GAMMA 1,GAMMA2) 

CALL S0LV4( GAMMA 1 , GAMMA2 ) 

UE(NX) = U(NP, 2) 

IF( ABS(DELV( 1) ) . GT. EPSL) GO TO 70 
C 

C CHECK ON LAMINAR FLOW SEPARATION. IF SEPARATION OCCURS, ASSIGN BEGIN 
C OF TRANSITION TO THAT POINT AND RECOMPUTE THE CURRENT STATION NX 
C 

IF(V(1,2). GT. 0. 0 .OR. ITR. NE.3) GOTO 110 
CALL TRNS( I CODE ) 

GOTO 25 
C 

C TURBULENT FLOW CALCULATION 
C 

100 CONTINUE 
CALL EDDY 

CALL COEF( GAMMA 1,GAMMA2) 

CALL S0LV4( GAMMA 1 , GAMMA 2 ) 

UE(NX) = U(NP, 2) 

VM = AMAX1(V(1,2) ,1. 0) 

IF(ABS(DELV(1)/VM) . GT. EPST) GO TO 70 
110 CONTINUE 
C 

C CHECK FOR B. L. GROWTH 
C 

IF(NP . GE. NP r ) GO TO 120 

IF( ABS(V(NP,2) ) . LT. 0.0005 .AND. ABS( 1. 0-U(NP-2,2)/U(NP,2)) 

+ . LT. 0.0035. OR. IGROW. GT. IGROWT) GOTO 120 

CALL FILLUPfn 
IGROW=IGROW+l 
IT ‘ = 1 

GO TO 70 
C 

120 CONTINUE 

CALL FILLUP( 2 ) 

CALL OUTPUT( 2 ) 

IF(NX. GE.NTR .OR. ITR. EQ. 0) GOTO 150 
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IF(NX. LT. 3 .OR. ITR. NE. 3) GOTO 150 
C 

C CALCULATE TRANSITION LOCATION USING MICHEL METHOD 
C 

CALL TRNS( I CODE) 

IF( ICODE. EQ. 0) GOTO 150 
C 

C TRANSITION OCCURS BASED ON MICHEL CRITERIOR AT STATION NX 

C RECALCULATE B. L. AT NX STATION ASSUMING THE FLOW IS TRANSITIONAL 

C 

IT =0 

I GROW = 1 
GOTO 70 

150 CONTINUE 

IF( . NOT. HOMOPY ) GO TO 154 
I F ( RT .GT. 0.9999) GO TO 154 
IRC = IRC + 1 
RT = RTSS(IRC) 

UEO(NX) =RT*UESAVE + ( 1. 0-RT)*UEREF 
GO TO 30 
154 CONTINUE 

IF (NX . LT. NXSTOP ) GO TO 20 
C 

C THE B. L. CALCULATION FOR THE CURRENT SWEEP IS COMPLETED. 

C CHECK FOR THE CONVERGENCE AND , IT NOT, MOVE TO THE NEXT 

C SWEEP. 

C 

c 

160 CONTINUE 

D(NXT) = GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,D(NXT-3) ,D(NXT-2) , 

+ D(NXT-l) ,X(NXT) ) 

DLS(NXT)= GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,DLS(NXT-3) ,DLS(NXT-2) , 
+ DLS(NXT-l) ,X(NXT) ) 

UE(NXT) = GRANG(X(NXT-3),X(NXT-2),X(NXT-l),UE(NXT-3), 

+ UE ( NXT - 2 ) , UE ( NXT - 1 ) , X ( NXT ) ) 

DO 165 I = 1 , NXSTOP 
165 CFS(I) = CF( I ) 

CALL AME AN( NS, NXSTOP, X,CF,1) 

C CALL AME AN( NS, NXSTOP, X,VW,1) 

C CALL HEADER( TITLE , SURFID, ISTRP ) 

IF ( ICYCLE .LT. ICYTL-1 .AND. IP . LT. 0)G0 TO 170 
WRITE (6, 250 ) ISWP 

WRITE (6, 262 ) 1 ,XC( 1) ,X( 1) ,CF( 1) ,DLS( 1) ,THT( 1) ,UE( 1) , 

+ UE0( 1) ,0. 0,GMTRS( 1) , ITP( 1) ,NPSTR( 1) 

WRITE (6, 264 ) (M,XC(M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 

+ UEO(M) ,DLS(M)/THT(M) ,GMTRS(M) , ITP(M) ,NPSTR(M) ,M=2, NXSTOP) 

IF ( (ICYCLE. EQ. ICYTL). AND. (IP. EQ. -2). AND. (NVP(ISF).NE. 0)) THEN 
WRITE( 12,800) NS+ 1 , NTR , XCTR 

WRITE ( 12,810) (XC(M) ,X(M) ,UE(M) ,CF(M) ,GMTRS(M) , ISG(M) , 

2 M=l, NXSTOP) 

800 FORMAT( 215, F10. 6) 

810 FORMAT(5E15. 5,15) 

END IF 

170 CONTINUE 
C 

DMAX = D(l) 
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180 

C 

C 

C 

C 



190 

195 

200 

205 



C 

C 

C 



DDMAX = ABS(D(1) - DB(1)) 

DO ISO I = 2 , NXT 
DM AX = AMAX1( DMAX,D(I) ) 

DD = ABS( D( I ) - DB( I ) ) 
DDMAX = AMAX1( DDMAX, DD ) 
CONTINUE 

IF ( ABS( DDMAX / DMAX ) . LE. 



UPDATE D FOR THE NEXT SWEEP 



0.0050 ) RETURN 



210 

C ■ 

250 

260 

262 



264 

270 



IF (ISWP . GT. 1) GO TO 195 
DO 190 I = NS , NXT 

DC I) = D( I )*( 1. 0+0MEGA*(UE( I) /UEO( I ) - 1. 0) ) 

GO TO 205 

IF (ISWP .EQ. 2) GOTO 205 
DO 200 I = NS , NXT 

D( I) = D( I) * (1. 0+0MEGA*(UE(I)/UEB(I)-l. 0)) 

IF( ISWP . GE. ISWPT ) RETURN 
NX = NS 

NP = NPSTR(NX) 

INDEX = 2 

DO 210 1= 1 , NXT 

DB(I) = D( I) 

UEB(I) = UE( I ) 

CONTINUE 
GOTO 10 

FORMAT( 1H0 , ’ *-v SUMMARY OF INVERSE BOUNDARY LAYER SOLUTIONS. **’ , 
+ 1H0 , 4X, ' ISWP =' , 13/) 

FORMAT ( 1H0 , 4X , 2HNX , 5X , 3HX/C , 9X , 1HX , 9X , 2HCF , 8X , 

+ 3HDLS , 8X , 3HTHT , 9X , 2HUE , 8X,3HUE0, 10X, 1HD ,9X,2HDB,3X, 

+ 4HGMTR , 4X , 2HIT , IX , 2HNP/ ( 1H , 3X, 13 ,F10. 5 , 8E11. 4 ,F8. 4, 213) ) 

FORMAT ( 1H0 , 4X , 2HNX , 6X , 3HX/C , 1 IX , 1HX , 10X , 2HCF , 9X , 

+ 3HDLS , 9X , 3HTHT , 10X , 2 HUE , 9X , 3HUE0 , 1 IX , 1HH , 8X , 

+ 4HGMTR , 4X , 2HIT , IX , 2HNP/ ( 1H ,3X, 13, 9E12. 4,213)) 

FORMAT ( 1H , 3X, 13, 9E12. 4,213) 

FORMATC 1H0 , ' ** ITERATIONS EXCEEDED ITMAX AT NX =’ , 15 / 

+ 1H0 , 1 ** CALCULATIONS STOP. **’) 

END 

DATA SET KCBCOUTPUT AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCOUTPUT AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCOUTPUT AT LEVEL 002 AS OF 02/22/84 

SUBROUTINE OUTPUT( INDEX) 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON/BLIN/ TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR,XTR, ISTRP, ICYCLE , ICYTL,XCTRS( 2) ,TRFIND( 2) 

COMMON /BLC 1/ F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON /BLC7/ C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
C0MM0N/EDDY1/ RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) 

+ , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON /GRD / ETA( 101) ,DETA( 101) ,A( 101) 

COMMON /SMRY/ VW( 100) ,ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) ,THT( 100) 
+ NPSTR(IOO) 

COMMON /I SURF/ ISF 



INT22260 
INT22270 
INT22280 
INT22290 
INT22300 
INT22310 
INT22320 
INT22330 
INT22340 
INT22350 
INT22360 
INT22370 
INT22380 
INT22390 
INT22400 
INT22410 
INT22420 
INT22430 
INT22440 
INT22450 
INT22460 
INT22470 
INT22480 
INT22490 
INT22500 
INT225 10 
INT22520 
INT22530 
/INT22540 
INT22550 
INT22560 
INT22570 
INT22580 
INT22590 
INT22600 
INT22610 
INT22620 
INT22630 
INT22640 
INT22650 
INT22660 
INT22670 
INT22680 
INT22690 
INT22700 
INT227 10 
INT22720 
INT22730 
INT22740 
INT22750 
INT22760 
INT22770 
INT22780 
, INT22790 
INT22800 
INT22810 
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c 

p 


COMMON/PLOT/NVP( 2) ,NXVP( 20 , 2) , ICC 


INT22820 

INT22830 

TVT99RAn 


o 






c 




INT22850 




ITP(NX) = IT 


INT22860 




NP STR ( NX ) =NP 


INT2287 0 




IF(NX. GT. 1) GOTO 5 


INT22880 




DLS(NX)= 0. 0 


INT22890 




VW(NX) = 0. 0 


INT22900 




D( NX) = 0. 0 


INT229 10 




THT( NX)= 0.0 


INT22920 




CF(NX) = 0. 0 


INT22930 




W(NX) = 0. 0 


INT22940 




GOTO 150 


INT22950 


5 


GOTO (10,100,200), INDEX 


INT22960 


c 




INT22970 


c 


CALCULATE B. L. PARAMETERS FOR TRANSFORMED COORDINATES 


INT22980 


10 


CONTINUE 


INT22990 




CF(NX) =2.0 * V( 1 , 2) * B( 1 ,2)/SQRX 


INT23000 




W(NX) = UE(NX)*SQRT(UE(NX)/X(NX))*V(1,2) 


INT23010 




DLS(NX)= X(NX)/SQRX * (ETA(NP) -F(NP, 2) ) 


INT23020 




D(NX) = UE(NX) * DLS(NX) * SQRT(RL) 


INT23030 




U1 = U(l,2) * (1. 0 -U( 1,2)) 


INT23040 




SUM = 0. 0 


INT23050 




DO 20 J=2 ,NP 


INT23060 




U2 = U(J,2) * ( 1. 0 -U( J , 2) ) 


INT23070 




SUM = SUM + A(J) * (U1 + U2) 


INT23080 




U1 = U2 


INT23090 


20 


CONTINUE 


INT23100 




THT(NX)= X(NX)/SQRX * SUM 


INT231 10 




GOTO 150 


INT23120 


C 




INT23130 


C 


CALCULATE B. L. PARAMETERS FOR SEMI-TRANSF COORDINATES 


INT23140 


100 


CONTINUE 


INT23150 




SQXC = SQRT(X(NX) ) 


INT23160 




SQRL = SQRT(RL) 


INT23170 




CF(NX) =2.0 * V( 1 , 2) * B ( 1 , 2 ) / ( SQXC*SQRL*W ( NP , 2 ) **2 ) 


INT23180 




VW(NX) = V( 1,2) / SQXC 


INT23190 




UE(NX) = U(NP, 2) 


INT23200 




RX = RL * UE(NX) * X(NX) 


INT23210 




DLS(NX) = (ETA(NP) -F(NP,2)/U(NP, 2) )/SQRL*SQXC 


INT23220 




SUM =0.0 


INT23230 




U1 = U(1,2)/U(NP,2)*(1. 0 -U( 1 , 2)/U( NP, 2) ) 


INT23240 




DO 120 J=2 ,NP 


INT23250 




U2 = U( J,2)/U(NP,2)*( 1. 0 -U( J, 2)/U(NP, 2) ) 


INT23260 




SUM = SUM + A( J) * (U1 + U2) 


INT23270 




U1 = U2 


INT23280 


120 


CONTINUE 


INT23290 




THT(NX) = SUM /SQRL * SQXC 


INT23300 




D(NX) = ( U( NP , 2 ) *ETA( NP)-F(NP,2)) * SQXC 


INT23310 


150 


IF (NX . GE. NXT) GO TO 160 


INT23320 




IF ( IEDY .EQ. 0 .OR. NX . LE. NTR+2 ) GO TO 160 


INT23330 


C 




INT23340 


c : 


MODIFY ALFA USING SIMPSON'S ARGUMENTS 


INT23350 


c 




INT23360 



105 





CALL SMPSON 


INT23370 


160 


DO 175 J=1 , NPT 


INT23380 




F( J, 1) = F(J,2) 


INT23390 




U( J, 1 ) = U(J 2) 


INT23400 




V ( J , 1 ) = V(J,2) 


INT23410 




W(J,1) = W(J,2) 


INT23420 




B( J, 1) = B ( J . 2 1 


INT23430 


175 


CONTINUE 


INT23440 




IF ((IP. LE. 0). AND. ((IP.NE. -2). OR. (ICYCLE.LT. ICYTL))) RETURN 


INT23450 


C 




INT23460 


C 


PRINT OUT VELOCITY PROFILES 


INT23470 


200 


IF (NX. EQ. 1) GOTO 210 


INT23480 ) 




IF (NX.LE.NS) THEN 


INT23490 




FAC1 = SQRT(X(NX) /RL/UE(NX) ) 


INT23500 




FAC2 = 1.0 


INT235 10 




ELSE 


INT23520 




FAC1 = SQRT(X(NX)/RL) 


INT23530 




FAC2 = 1. 0/UE(NX) 


INT23540 




ENDIF 


INT23550 




NPM1 = NP -1 


INT23560 




WRITE(6 ,4001) NX,X(NX) 


INT23570 




WRITE( 6 ,4000) 


INT23580 




WRITE(6 ,4100) ( J,ETA( J) ,F(J,2),U(J,2),V(J,2),W(J,2),B(J,2), 


INT23590 




+ ETA( J) *FAC 1 ,U( J, 2 )*FAC2 , J=1 , NPM1 , 3 ) 


INT23600 



c 

210 



WRITE (6, 4 100) NP,ETA(NP),F(NP,2) ,U(NP,2) ,V(NP,2) ,W(NP,2) ,B(NP,2) , 
+ ETA(NP)*FAC1,U(NP,2)*FAC2 



IF (IP.NE. -2) RETURN 
IF ( (NXVP( ICC , ISF) . NE. NX) . OR. ( ICC. GT. NVP( ISF) ) ) RETURN 
WRITE( 12,4200) NP 
WRITE( 12,4300) (ETA( J) , J=1,NP) 

WRITE( 12,4300) (U( J,2) , J=1,NP) 

ICC = ICC+1 
RETURN 

4001 FORMAT(/1HO, 'NX =' ,15, ' S/C=’,F10.5) 

4000 FORMAT ( 1H0 , 2H J,9X, 3HETA , 15X , 1HF , 13X , 1HU , 13X , 1HV , 13X , 1HW , 1 3X , 1HB , 
+ 13X, 3HY/C , 10X, 4HU/UE) 

4100 FORMAT ( 1H , 13 ,E14. 5 , 2X,5E14. 5 , 2X, 2E14. 5) 

4200 FORMAT( 15 ) 

4300 FORMAT( 8F10. 6) 

END 

DATA SET KCBCSMFIT AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCSMFIT AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCSMFIT AT LEVEL 001 AS OF 08/15/83 

SUBROUTINE SMFIT( NS , ND , X , Q , D , KS ) 



C 

C 

C 

C 

C 

C 

C 



THIS SUBROUTINE SMOOTHES DATA, Q, USING FIVE -POINT FORMULA. 



NS 

X 

Q 

KS 



BEGINNING POINT? ND : END POINT 
INDEDENPENT COORDINATE? D : WORKING STORAGE 
VARIABLE TO BE SMOOTHED 
NO OF SMOOTHING 



DIMENSION X( 101) ,Q( 101) ,D( 101) 



SMT5(Q1,Q2,Q3,Q4,Q5 ) = 0. 0625*( 10. 0*Q3+4. 0*( Q2+Q4) -Q1-Q5 ) 

SMT3(Q1,Q2,Q3,X1,X2,X3) = 0. 5*(Q2+(Q1*ABS(X3-X2)+Q3*ABS(X2-X1)) 



INT23610 
INT23620 
INT23630 
INT23640 
INT23650 
INT23660 
INT23670 
INT23680 
INT23690 
INT23700 
INT237 10 
INT23720 
INT23730 
INT23740 
INT23750 
INT23760 
INT237 70 
INT23780 
INT23790 
INT23800 
INT23810 
INT23820 
INT23830 
INT23840 
INT23850 
INT23860 
INT23870 
INT23880 
INT23890 
INT23900 
INT239 10 
INT23920 
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20 

40 

100 

C 

200 



250 

300 

C 



r /ABS(X3-X1) ) 

IF(KS. LE. 0) RETURN 

NSP1 = NS+1 

NSP2 = NS+2 

NDM1 = ND-1 

NDM2 = ND-2 

NDIF = ND-NS+1 

IF ( NDIF .LT. 3 ) RETURN 

IF ( NDIF . LT. 5 ) GO TO 200 

DO 100 K=1 ,KS 

D(NS+1)= SMT3f Q(NS) ,Q(NS+1) ,Q(NS+2) ,X(NS) ,X(NS+1) ,X(NS+2) ) 
D(ND-1)= SMT3' 0(ND-2) ,Q(ND-1) ,Q(ND) ,X(ND-2) ,X(ND-1) ,X(ND)) 
DO 20 I=NSP2,NDM2 

D( I) = SMT5(Q( 1-2) ,Q( 1-1) ,Q( I) ,Q( 1+1) , Q( 1+2) ) 

CONTINUE 

DO 40 I=NSP i , NDM1 
Q( I) = DC I) 

CONTINUE 

RETURN 



DO 300 K = 1 , KS 
DO 220 I = NSP1 ,NDM1 

D( I) = SMT3(Q(I-1),Q(I),Q(I+1),X(I-1),X(I),X(I+1)) 

220 CONTINUE 

DO 250 I = NSP1 ,NDM1 
Q( I) = DC I) 

CONTINUE 
CONTINUE 

RETURN 

END 

C DATA SET KCBCS0LV3 AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCS0LV3 AT LEVEL 001 AS OF 08/24/84 

C DATA SET KCBCS0LV3 AT LEVEL 005 AS OF 02/21/84 

SUBROUTINE S0LV3 

COMMON /BLCO/ NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 

COMMON /BLC 1/ F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON /BLC2/ DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON /BLC6/ Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
+ S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON /GRD / ETA( 101) ,DETA( 101) , A( 101) 

COMMON /BLCB/ Al 1( 101) ,A12( 101) , A13( 101) , A14( 101) , 

+ A21( 101) , A22( 101) , A23( 101) ,A24(101) 

C 

A 1 1 ( 1)= 1. 0 
A12( 1)= 0. 0 
A13( 1 )= 0. 0 
A21( 1)= 0. 0 
A22( 1)= 1. 0 
A23( 1)= 0. 0 
Gll=- 1. 0 
G12=-A( 2) 

G13= 0. 0 



INT23930 
INT23940 
INT23950 
INT23960 
INT23970 
INT23980 
INT23990 
INT24000 
INT240 10 
INT24020 
INT24030 
INT24040 
INT24050 
INT24060 
INT24070 
INT24080 
INT24090 
INT24100 
INT241 10 
INT24120 
1NT24130 
INT24140 
INT24150 
INT24160 
INT24170 
INT24180 
INT24190 
INT24200 
INT24210 
INT24220 
INT24230 
INT24240 
INT24250 
INT24260 
INT24270 
INT24280 
INT24290 
INT24300 
INT243 10 
INT24320 
INT24330 
INT24340 
101) INT24350 
INT24360 
INT24370 
INT24380 
- - INT24390 
INT24400 
INT24410 
INT24420 
INT24430 
INT24440 
INT24450 
INT24460 
INT24470 
INT24480 
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c 

c 

c 



G21= S4( 2) 

G23=-S2(2)/A(2) 

G22= G23+S6( 2 ) 

FORWARD SWEEP 

DO 500 J=2 ,NP 

IF( J . EQ. 2) GO TO 100 

DEN = (A13(J-1)*A21(J-1)-A23(J-1)*A11(J-1)-A(J)* 

+ (A12( J-1)*A21(J-1)-A22( J-1)*A11(J-1))) 

DENI = A22( J- 1)*A( J) -A23( J-l ) 

G1 1= (A23( J-l)+A( J)*(A(J)*A21( J-l)-A22( J-1)))/DEN 
G12= -( A( J)*A( J)+G11*( A12( J-l )*A( J) -A13( J-l) ) )/DENl 
G13= (G11*A13(J-1)+G12*A23CJ-1))/AC J) 

G21= (S2(J)*A21(J-1)-S4(J)*A23(J-1)+A(J)*(S4(J)* 

+ A22( J-l)-S6( J)*A21( J-1)))/DEN 

G22= ( -S2( J)+S6( J)*A( J) -G21*( A( J)*A12( J-l) -A13( J-l) ) )/DENl 
G23= G21*A12( J-1)+G22*A22( J-1)-S6(J) 

100 A11(J)= 1. 0 

A12( J)=-A( J) -G13 
A 1 3 ( J ) = A(J)*G13 
S3( J) 

S5( J) -G23 
S1(J)+A(J)*G23 

R 1 ( J ) - ( G 1 1*R1 ( J - 1 ) +tj 1 2*R2 ( J - 1 ) +G 1 3*R3 ( J - 1 ) ) 

R2( J) G21*R1( J-1)+G22*R2( J-1)+G23*R3( J-l) ) 



500 
C 
C 
C 



A21( J)= 
A22( J)= 
A23( J)= 

RlCJ) = 

R2(J) = 
CONTINUE 



BACKWARD SWEEP 

DELUfNP) 

El 
E2 

DELV(NP) 



600 



DELF(NP) 

J 

J 

E3 

DEN2 = 

h 

DELV(J) 

h 

DELU(J) 
DELF( J) 
IF( J . GT 



NP 

J-l 



R3(NP) 

Rl(NP) -A12(NP)*DELU(NP) 

R 2 ( NP ) - A2 2 ( NP ) *DE LU ( NP ) 

(E2*A11(NP)-E1*A21(NP))/(A23(NP)*A11(NP)-A13(NP)* 
A21(NP) ) 

(E1-A13(NP)*DELV(NP))/A11(NP) 



INT24490 
INT24500 
INT245 10 
INT24520 
INT24530 
INT24540 
INT24550 
INT24560 
INT24570 
INT24580 
INT24590 
INT24600 
INT24610 
INT24620 
INT24630 
INT24640 
INT24650 
INT24660 
INT24670 
INT24680 
INT24690 
INT24700 
INT247 10 
INT24720 
INT24730 
INT24740 
INT24750 
INT24760 
INT24770 
INT24780 
INT24790 
INT24800 
INT24810 
INT24820 
INT24830 
INT24840 
INT24850 
INT24860 
INT24870 
INT24880 
INT24S90 



R3(J)-DELU(J+1)+A(J+1)*DELV(J+1) 

A21(J)*A12(J)*A(J+1)-A21(J)*A13(J)-A(J+1)*A22(J)*A11(J)+ 

A23( J)*A11( J) 

= (All( J)*(R2(J)+E3*A22(J))-A21(J)*R1(J)-E3*A21( J)*A12( J) INT24900 
)/DEN2 INT249 10 

=-A( J+1)*DELV( J) -E3 

= (Rl( J)-A12( J)*DELU( J)-A13(J)*DELV( J))/A11(J) 

1) GO TO 600 



700 



DO 700 J=1 ,NP 
F(J,2)= F( J,2)+DELF( J) 

U( J, 2)= U( J , 2)+DELU( J) 

V(J,2)= V( J,2)+DELV( J) 

CONTINUE 
U(l,2)= 0. 0 

CALL EDGCHK(NP,ETA,F(1,2),U(1,2),V(1,2)) 
RETURN 



INT24920 
INT24930 
INT24940 
INT24950 
INT24960 
INT249 70 
INT24980 
INT24990 
INT25000 
INT25010 
INT25020 
INT25030 
INT25040 
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o o o o o o 



10 

c 

c 



END 

DATA SET KCBCS0LV4 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCS0LV4 AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCS0LV4 AT LEVEL 001 AS OF 02/21/84 

SUBROUTINE S0LV4( GAMMA 1 ,GAMMA2) 

COMMON /BLCO/ NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 

COMMON /BLC1/ F( 101 ,2) ,U( 101 , 2) ,V( 101 , 2) ,W( 101 , 2) ,B( 101 ,2) 

COMMON /BLC2/ DELF( 101) ,DEL T J( 101) ,DELV( 101) ,DELW( 101) 

COMMON /BLC6/ Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 



INT25050 
INT25060 
INT25070 
INT25080 
INT25090 
INT25100 
INT25110 
INT25120 
INT25130 
101) INT25 140 



COMMON 


/GRD / ETA( 101) ,DETA( 101), A( 101) 


INT25150 


COMMON 


/BLCB/ All( 101) , A12( 101) ,A13(101) ,A14(101) , 


INT25160 




A21(101),A22(101),A23(101),A24(101) 


INT25170 
INT25180 
- - - INT25190 
INT25200 


A 1 1 ( 1 ) 


= 1.0 


INT252 10 


A12( 1) 


= 0.0 


INT25220 


A13( 1) 


= 0.0 


INT25230 


A14( 1) 


= 0.0 


INT25240 


A21( 1) 


= 0. 0 


INT25250 


A22( 1) 


= 1. 0 


INT25260 


A23( 1) 


= 0. 0 


INT25270 


A24( 1) 


= 0.0 


INT25280 


DO 10 J 


= 2 ,NP 


INT25290 


AAl 


= Al3(J-l) -A( J)' V A12( J-l) 


INT25300 


AA2 


= A23( J-l) -A( J) *A22( J- 1) 


INT25310 


AA3 


= S2(J)-A(J)*S6(J) 


INT25320 


DET 


= AA2*A11( J-l) -AA1*A21( J-l) 


INT25330 


AJS 


= A(J)**2 


INT25340 


Gil 


-( AA2+A2 1( J- 1) ,V AJS) /DET 


INT25350 


G12 


( A 1 1 ( J - 1 ) * A JS+A A 1 ) /DET 


INT25360 


G13 


A12(J-1)*G11+A22(J-1)*G12+A(J) 


INT25370 


G14 


A14(J-1)*G11+A24(J-1)*G12 


INT25380 


G21 


( S4( J)*AA2 -A2 1( J- 1)*AA3) /DET 


INT25390 


G22 


(All(J-l) - V AA3 - S4( J ) *AA 1 ) /DET 


INT25400 


G23 


A12(J - I)*G21+A22(J-1)*G22-S6(J) 


INT25410 


G24 


A14( J- 1)*G21+A24( J-1)*G22-S8( J) 


INT25420 


All(J) 


= 1. 0 


INT25430 


A12( J) 


= ') -G13 


INT25440 


A13( J) 


= A(J)*G13 


INT25450 


A14( J) 


= -G14 


INT25460 


A21( J) 


= S3( J) 


INT25470 


A22( J) 


= S5( J) -G23 


INT25480 


A23( J) 


= S1(J)+A(J)*G23 


INT25490 


A24( J) 


= S7( J) -G24 


INT25500 


R1(J) 


= Rl( J) -G1 1*R1( J- 1) -G12*R2( J- 1) -R3( J- 1)*G13 


INT255 10 




-G14*R4( J- 1) 


INT25520 


R2(J) 


= R2( J) -G2 1*R1( J- 1) -G22*R2( J- 1) -R3( J-1)*G23 


INT25530 




-G24*R4( J- 1) 


INT25540 


CONTINUE 


INT25550 






INT25560 


iCKWARD 


SWEEP 


INT25570 


J 


= NP 


INT255S0 


G1 


= GAMMA 1/GAMMA2 


INT25590 


R3( J) 


= R3( J)/GAMMA2 


INT25600 
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n o o n n o 



20 



30 



C 



C 

100 

110 

C 



R1(J) 

R2(J) 

Cl 
C2 
DET 

DELF(J) 

DELV(J) 

DELW(J) 

DELU( J) 

J 

CC1 
CC2 
CC3 
CC4 
CC5 
CC6 
DENO 
DELF( J) 

DELV(J) 

DELW( J) 

DELU(J) 

IF( J . GE. 

DO 30 J = 1 ,NP 
F( J , 2 ) = F(J,2)+DELF(J) 
U( J, 2)+DELU( J) 
V( J,2)+DELV(J) 
W(J,2)+DELW(J) 



= Rl( J)-A12( J)*(R4( J)+R3( J))-A14(J)*R3( J) 
= R2(J)-A22(J)*(R4( J)+R3( J))-A24( J)*R3( J) 
- A11(J)-G1*(A12(J)+A14(J)) 

= A21(J)-G1*(A22(J)+A24(J)) 

= C1*A23(J)-C2*A13(J) 

= (R1(J)*A23(J)-R2(J)*A13(J))/DET 
= (C1*R2(J)-C2*R1(J))/DET 
= R3(J)-G1*DELF(J) 

= R4(J)+DELV(J) 

= J-l 

= DELU( J+l) -R3( J) -A( J+1)*DELV( J+l) 

= DELW( J+l) -R4( J) 

= A13(J)-A(J+1)*A12(J) 

■ Rl( J)-A12( J)*CC1-A14(J)*CC2 
= A23(J)-A(J+1)*A22(J) 

= R2( J)-A22( J)*CC1-A24(J)*CC2 
= A11(J)*CC5-A21( J)*CC3 
= ( CC4*CC5 -CC3*CC6 ) /DENO 
= ( A 1 1 ( J ) *CC6 - A2 1 ( J ) *CC4 ) /DENO 
= CC2 

= CC1-A(J+1)*DELV(J) 

2) GO TO 20 



U(J,2) = 
V(J,2) = 
W(J,2) = 
CONTINUE 
U(l,2) = 



0 . 0 



CALL EDGCHK(NP,ETA,F( 1,2) ,U( 1,2) ,V(1,2)) 
RETURN 



END 

DATA SET KCBCTRNS AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCTRNS AT LEVEL 001 AS OF 08/24/84 

DATA SET KCBCTRNS AT LEVEL 005 AS OF 03/13/84 

SUBROUTINE TRNS(ICODE) 

CALCULATE TRANSITION LOCATION USING MICHEL CRITERION 

COMMON /BLCO/ NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

COMMON /BLC1/ F( 101,2) ,U(101,2) ,V(101 ,2) ,W(101 ,2) ,B( 101,2) 

COMMON/BLIN/ TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+ XCTR , XTR , I STRP , I C YC LE , I C YTL , XCTRS ( 2 ) ,TRFIND( 2) 

COMMON/EDDY 1 / RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) 

+ ,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON /GRD / ETA( 101) ,DETA( 101) , A( 101) 

COMMON /GTY / X( 101) ,UE( 100) ,P1( 100) , P2( 100) ,CEL,CELH 
COMMON /TRN/ PGAMTR , OMEGA , RTHETB , RTRANB 



F0RMAT(/3X. 'ELGIN OF TRANSITION HAS BEEN DETECTED BY MICHEL"S ' 
+ ' CRITERION: X/C =' ,F8. 4,4X, 1 S/C =' ,F8. 4 ,4X, ' NTR =’ ,13/) 
F0RMAT(/3X, ’BEGIN OF TRANSITION IS ASSUMED AT THE POINT OF 
+' LAMINAR SEPARATION: X/C =' ,F8. 4,4X, ' S/C =' ,F8. 4,4X, ’ NTR =' ,13/) 



I CODE = 0 
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ISEP = 1 

c 

IF( V( 1,2). LT. 0. 0) THEN 

C *** TRANSITION PROCESS HAS BEGUN DUE TO LAMINAR SEPARATION *** 

FAC = V(1,1)/(V(1,1)-V(1,2)) 

GOTO 20 
END IF 

*** CHECK MICHEL'S TRANSITION CRITERION *** 

ISEP = 0 

SUM = 0. 0 

FI = U( 1 , 2)/U(NP , 2)*( 1. 0-U( 1 , 2)/U(NP, 2) ) 

DO 10 J=2 ,NP 

F2 = U( J, 2)/U(NP , 2)*( 1. 0-U(J,2)/U(NP,2)) 

SUM = SUM + (FI + F2 ) *A(J) 

FI = F2 

10 CONTINUE 

CONV = SQRT(RL/X(NX) ) 

IF(NX. LE.NS) CONV = SQRT(RX)/X(NX) 

THETA = SUM / CONV 
RTHETA = RL * UE(NX) * THETA 
RTRAN = 1.174 * ( 1. 0+22400. 0/RX) * RX**0. 46 
IF( RTHETA. LT. RTRAN) THEN 
RTHETB = RTHETA 
RTRANB = RTRAN 
RETURN 
END IF 
C 

C *** TRANSITION PROCESS HAS BEGUN BECAUSE OF MICHEL'S CRITERION *** 
FAC = (RTHETB -RTRANB )/ (RTRAN - RTRANB -RTHETA+RTHETB) 

C 

C *** COMPUTE EXACT LOCATION OF TRANSITION BEGIN *** 

20 NTR = NX-1 
NTR1 = NTR + 1 

XCTR = XC(NX-l) + FAC*(XC(NX) -XC(NX-l) ) 

XTR = X(NX-l) + FAC*(X(NX) -X(NX- 1) ) 

UETR = UE(NX-l) + FAC*(UE(NX) -UE(NX-l) ) 

IF ( ISEP . EQ. 0 ) WRITE (6,100) XCTR, XTR, NTR 
IF ( ISEP .EQ. 1 ) WRITE (6,110) XCTR, XTR, NTR 
ICODE = 1 
C 

C *** CALCULATE I NTERM ITTENCY DISTRIBUTION *** 

RXNTR = XTR * UETR * RL 
GGFT = RL**2/PGAMTR/ RXNTR** 1. 34*UETR**3 
DO 30 I=NTR1 ,NXT 
ALFAS(I) = 0. 0168 
GMTRS( I )= 1. 0 
30 CONTINUE 

ALFAS(NTR) = 0.0168 
UEINTG = 0. 0 
U1 = 0.5/UETR 
XI = XTR 

DO 40 I=NTR1,NXT 
U2 = 0. 5/UE( I) 

X2 = X(I) 

UEINTG = UEINTG+(U1+U2)*(X2-X1) 
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Ul = U2 




INT26730 




XI = X2 




INT26740 




GG = GGFT*UEINTG*(X(I)-XTR) 




INT26750 




IF( GG . GT. 10.0) GOTO 50 




INT267 60 




GMTRS(I) = l.O-EXP(-GG) 




INT26770 




40 CONTINUE 




INT26780 


c 






INT26790 


C 


•** RESET FINITE DIFFERENCE CALCULATIONS *** 




INT26800 




50 DO 60 J=1 ,NPT 




INT268 10 




F(J,2) = F(J,1) 




INT26820 




U(J,2) = U(J,1) 




INT26830 




V( J ,2) = V(J,1) 




INT26840 




B(J,2) = B(J,1) 




INT26850 




W(J,2) = W ( J , 1 ) 




INT26860 


60 


CONTINUE 




INT26870 




RETURN 




INT26880 




END 




INT26890 


C 






INT26900 


c 


SUBROUTINE EDGCHK(NP, ETA, F, U, V) 




INT269 10 
INT26920 


c 


DIMENSION ETA(lOl), F(101), U(101), V(101) 




INT26930 

INT26940 


c 


JS = NP - 3 




INT26950 

INT26960 




NPM1 = NP - 1 




INT269 70 




DO 10 J=JS , NPM1 




INT26980 




JJ = J 




INT26990 




IF(U(J). GE.U(NP) .OR. V(J).LT. 0.0) GOTO 20 




INT27000 


10 


CONTINUE 




INT27010 




RETURN 




INT27020 


20 


JS = JJ - 1 




INT27030 




IF( JS. GT. (NP-2) ) JS = NP-2 




INT27040 




CALL AMEAN( JS, NP, ETA, U, 1) 




INT27050 




CALL AMEAN(JS, NP, ETA, F, 1) 




INT27060 




DETAP = ETA(JS) -ETA(JS-l) 




INT27070 




VJP = (U(JS)-U(JS-1))/DETAP 




INT27080 




DO 30 J=JS ,NPM1 




INT27090 




DETAM = ETA( J+1)-ETA( J) 




INT27100 




VJM = (U(J+1)-U(J)) /DETAM 




INT27110 




V(J) = ( VJM “DETAP + VJP*DETAM) / ( DETAP+DETAM) 




INT27120 




VJP = VJM 




INT27130 




DETAP = DETAM 




INT27140 


30 


CONTINUE 




INT27150 




V(NP) = -V(NP-l) +2.0 * VJP 




INT27160 


C 


RETURN 

/V /V 0\ 0\ 0\ /V 0\ 0\ 0\ 0\ 0\ 0\ 0\ 0\ 0\ / I /V 0% 0\ 0\ ✓ V 0\ *\ 0\ /> #V 0\ /V 0\ 0\ 0\ 0\ 0\ 0\ 0\ 0\ 4\ 0\ n /% 0\ 0S 0\ 
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C 


NOTES: (FOR CHANGING FROM THE ORIGINAL PROGRAM) 


>V 


INT27190 


C 




* 


INT27200 


C 


1. 'EDDY' HAS BEEN MODIFIED BY ADDING 'FINT'. 


* 


INT27210 


C 


2. SUBROUTINE ’EDGCHK' HAS BEEN ADDED. 


/V 


INT27220 


C 


3. GROWTH LIMIT HAS BEEN ADDED FOR 2 IN ' MAIN2' . 


* 


INT27230 


C 




0\ 


INT27240 


C 

n 


^ m3 m *3* fcU 3 m 


Vr 


INT27250 

INT27260 

INT27270 


0 


0\ 0\ 0\ 0\ 0\ 0\ 0\ 0% 0 % 0 % 0% 0\ 0\ 0\ 0\ 0\ 0\ 0% 0\ 0\ 0\ 0\ 0\ 0S 0\ 0\ 0\ 0\ 0\ 0\ 0\ 0\ 4\ 4S 0\ 0\ 4\ /» 0\ 0\ 0\ 4\ 0\ 0\ 0\ 0\ 0S 0* 0\ 0\ 0\ 4% 0\ 0\ 0\ 0 \ 

END 
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SUBROUTINE SMPSON INT27280 

C INT27290 

COMMON /BLCO/ NX ,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP INT27300 

COMMON /BLC1/ F( 101 ,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101 ,2) INT27310 

COMMON /BLC7/ C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI INT27320 

COMMON/ EDDY 1 / RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) INT27330 

+ ,ALFAS(100),FFS(100),RTS(100),IEDY,NXSPT INT27340 

COMMON /GTY / X( 101) ,UE( 100) , Pl( 100) ,P2( 100) ,CEL,CELH INT27350 

COMMON /GRD / ETA(lOl) ,DETA( 101) ,A( 101) INT27360 

DIMENSION CRD( 12) ,RTD( 12) INT27370 

DATA RTD/O. 00,0. 05,0. 12,0. 20,0. 30,0. 40,0. 50,0. 60,0. 70, INT27380 

+ 0.80,0.90,1.00/ INT27390 

DATA CRD/5. 00,4. 75,4. 35,3. 80,3. 25,2. 85,2. 58,2. 37,2. 25, INT27400 

+ 2.15,2.06,2.00/ INT27410 

C INT27420 

C INT27430 

C INT27440 

C STEP 1 CALCULATE (DU/DX)/(DU/DY) INT27450 

C IF(NX. LT.NXSPT) GOTO 10 INT27460 

C INT27470 

C IN THE SEPARATED REGION, ALFA SET TO BE CONSTANT 1NT27480 

C ALFAS(NX)= ALFASP INT27490 

C RETURN INT27500 

C( INT275 10 

C INT27520 

C 10 CONTINUE INT27530 

C IF(V( 1,2). GT. 0.0) GOTO 20 INT27540 

C INT27550 

C SEPARATION OCCURS. ALFA SET TO BE THE PREVIOUS ITERATED VALUE INT27560 

C ALFASP = ALFAS(NX) INT27570 

C NXSPT = NX INT27580 

C RETURN INT27590 

C( INT27600 

C MODIFY OUTER EDDY BASED ON SIMPSON SUGGESTION INT27610 

TM =0.0 INT27620 

JM = 1 INT27630 

DO 30 J=2 ,NP INT27640 

TS = (B(J,2)-1. 0)* V( J, 2) INT27650 

IF(TS.LT. TM) GOTO 30 INT27660 

TM = TS INT27670 

JM = J INT27680 

30 CONTINUE INT27690 

VN'XM = 0. 5*(V(JM,2)+V(JM,1)) INT27700 

IF (NX . LE. NS) GOTO 35 INT27710 

DUDX = (U( JM, 2) -U( JM, 1) ) / (X(NX) -X(NX-l) ) INT27720 

GO TO 38 INT27730 

35 DUDX = CEL*(U( JM, 2) -U( JM, 1) )+P2(NX)*U( JM, 2)+0. 5*ETA( JM)* INT27740 

+ VNXM*( P2 ( NX) - 1 . 0 ) INT2 7 7 5 0 

38 RU = RL INT27760 

IF(NX. LE. NS)RU = RL * UEO(NX) * X(NX) INT27770 

RL2 = SQRT(RU) INT27780 

RR = DUDX/VNXM/RL2 INT27790 

C INT27800 

C STEP 2 : CALCULATE (UU - VV)/UV INT27810 

VNXM = 0.5*(V(1,2)+V(1,1)) INT27820 
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60 

c 

c 

70 



RT = VNXM/TM 

PRINT' (3X, 215, 3F10. 3)' , NX , JM , VNXM , TM , RT 
IF (RT . LT. 0. 0) RT = 0. 0 
IF(RT. GT. 1. 0) GOTO 60 

CR = 6.0 /( 1. 0 + 2.0 * RT*( 2. 0 -RT)) 

CR = 2. 0 

DO 40 1=2,12 

IF(RT. LT. RTD(i)) GOTO 50 
40 CONTINUE 
GOTO 70 

50 CR = CRD( I -i)+(CRD( I) -CRD( 1-1) )*(RT-RTD( I - 1) )/(RTD( I ) -RTD( I ■ 
GOTO 70 

CR = ( 1. 0 + RT ) /RT 

STEP 3 : CALCULATE FF 
FR = CR * RR 
IF(FR . GT. 0. 35) FR = 0. 35 
IF (FR . LT. -0. 8) FR = -0. 8 
FFS(NX)= (FFS(NX) + (1.0 -FR))/ 2.0 

RTS(NX)= RT 

ALFAS(NX)= 0. 0168/FFS(NX)**2. 5 
RETURN 



C( 

c 



10 

20 

30 



40 



END 

SUBROUTINE XSPACE( NI , NRITE , XI I , XLLT , RAD , NL1 , NR1 ) 
DIMENSION XI I ( 200) ,T( 200) 

DATA PI/3. 14159265359879/ 

RAD = PI 

NLEFT=NI-1 -NRITE 
NR4=NRITE/2 

IF((NRITE/2*2) .NE. NRITE) NR4=(NRITE+l)/2 

NL1=NR4+1 

NL2=NR4+NLEFT 

NR1=NL2+1 

NR2=NI 

PI2=0. 5* PI 

RAD2=( PI -RAD)/2. 0+PI2 

RAD3=RAD2+RAD 

SRT =RAD2/FL0AT(NR4) 

SRT2=SRT 

IF( (NRITE/2*2) . NE. NRITE) SRT2=RAD2 /FLOAT ( NR4 - 1 ) 
SLT = RAD/FLO AT(NLEFT) 

DO 10 1=1, NR4 

XII(I)=0. 5*( 1. 0+COS(FLOAT(I-1)*SRT)) 

DO 20 I=NL1 ,NL2 

XI I ( I )=0. 5 +X LLT* C 0 S ( F LO AT ( I -NL1 )*SLT+RAD2 ) 

DO 30 I=NR1 ,NR2 

XII(I)=0. 5*( 1. 0+COS(FLOAT(I-NR1)*SRT2+RAD3)) 
NA=(NI+l)/2 

IF((NI/2*2) .EQ. NI) NA=NI/2+l 
FN1=FL0AT(NA-1) 

FN2=FN1 

IF((NI/2*2) .EQ. NI) FN2=FLOAT(NA-2) 

DO 40 1=1, NA 
T(I)=FL0AT(NA-I)/FN1 
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50 



C 

C 



C 



C 

C 

C 

C 

C 

C 



60 



62 

65 



68 



C 



CALL AMEAN( 1 , NA , T , XI I , 1 ) 

XBIF = XII(l) - XII( 2) 

IF(XDIF . LT. 0.004) THEN 
BO 45 1=2,5 

XII(I) = XII( 1-1) -XDIF*3. 0 
CONTINUE 

CALL AMEAN(2,NA,T,XII ,10) 

END IF 

DO 50 I=NA,NI 
XII(I) = XII(NI-I+1) 

RETURN 

END 

SUBROUTINE TRGRID (N1 , XO , Y0,NI,NRITE,XLLT,N10,RAD,ID,NXSS) 
THIS SUB. IS TO REGRID SPACING NEAR TRAILING -EDGE 

DIMENSION X0( 200) ,YO(200) ,XI(200) ,YI(200) ,D1(200) ,D2(200) ,D3(200) 
+ XOOf 200) , Y00( 200) ,XII(200) ,YII(200) ,WX(200) ,WY(200) , 

+ WXI' 200) ,WYI(200) ,T(200) 



N20 



= N10 
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IF ( (Nl/2*2) . F.O Nl) N20 = N'10+1 
IF ( ( NI - ( N.T ,'Zj *2 ) . NE. 0) N1I= (NI-l)/2+l 
IF( (,‘NI - ( NI / 2 ) *2 ) .EQ. 0) NlI=NI/2 
N2I = Nil 

IF( (NI/2*2) .EQ. NI) N2I = N1I+1 

CALL XSPACE(NI ,NRITE ,XI ,XLLT,RAD,NL1 ,NR1) 
PRINT *, 1 NRITE=' ,NRITE , ' XLLT= ' ,XLLT 
WRITE (6, 290) 

WRITE (6, 300) (XO(I) ,1=1, Nl) 

WRITE (6, 298) 

WRITE (6, 300) (YO(I) ,I=1»N1) 

I F ( ID .EQ. 2) THEN 
DO 60 I=NL1 ,NXSS 
YI( I)=YO( I) 

XI(I)=XO(I) 

NXST=NXSS+8 

XM1=(XI(NXST) -XI(NXSS) )/8. 0 
XM2=(XI(NR1) -XI(NXSS) )/(NRl-NXSS) 

DO 62 I=NXSS ,NR1- 1 
XI(I)=(X0(I)+XI(I))/2. 0 
DO 65 I=NXSS ,NXST 
T( I ) =FLOAT( I -NXSS ) *XM1+XI ( NXSS ) 

CALL AMEAN(NXSS ,NXST,T,XI ,8) 

DO 68 I=NXSS,NR1 

T( I ) =FLOAT( I -NXSS )*XM2+XI ( NXSS ) 

CALL AMEAN(NXSS ,NR1 ,T,XI ,12) 

N10=NL1 
N1I=N10 
N20=N1+1-NXSS 
N2I=NI -NXSS+1 
ELSE 

NXSS=N2I 
END IF 

FOR LOWER SURFACE 
DO 5 I = 1 , N10 
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II = N10 - I + 1 

WX(II) = XO(I) 

WY (II) = YO( I) 

5 CONTINUE 

DO 7 I = 1 , Nil 

II = Nil - I + 1 

WXI(II) = XI(I) 

7 CONTINUE 

CALL DIFF3(N10,WX,WY,D1 ,D2 ,D3 ,0) 
CALL INTRP3(N10,WX,WY,D1,D2,D3,N1I, 
DO 9 I = 1 , Nil 

II = Nil - I + 1 

YI(II) = WYI(I) 

9 CONTINUE 
C 

C FOR UPPER SURFACE 

DO 10 I = 1 , N20 
II = N1 - N20 + I 
XOO(I) = XO(II) 

YOO(I) = YO(II) 

10 CONTINUE 

DO 20 I = 1 , N2I 
II = NI - N2I + I 
WXI(I) = XI(II) 

20 CONTINUE 

CALL DIFF3(N20,X00, YOO,Dl ,D2 ,D3, 0) 
CALL INTRP3 ( N20 , XOO , YOO , D1 , D2 , D3 , : 
DO 25 I = 1 , N2I 
II = N2I - I -• 1 
XII(II) = WXI(I) 

YII(II) = WYI(I) 

25 CONTINUE 

C 

C COMBINE TWO SURFACES INTO ONE CIRCLE 
C 

C NN = Nl - N10 - N20 
C DO 30 I = 1 , NN 

C II = Nil + I 

C 111= N10 + I 

C XI(II) = XO(III) 

C YI(II) = YO(III) 

C30 CONTINUE 

DO 40 I = 1 , N2I 

11 = NXSS + I -1 

12 = N2I - I + 1 
XI(I1) = XII( 12) 

Y I ( II) = YII( 12) 

40 CONTINUE 

XI(1) = XO(l) 

XI(I1)= XO(Nl) 

YI( 1) = YO( 1) 

YI( Il)= YO(Nl) 

C 

Nl = II 

DO 50 I = 1 , Nl 
XO(I) = XI(I) 
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YOCI) = YI( I) 


INT295 10 


50 


CONTINUE 


INT29520 


C 




INT29530 


C 


WRITE (6 , 295) 


INT29540 


C 


WRITE (6 , 300) (XO(I) , 1=1, Nl) 


INT29550 


C 


WRITE (6 , 298) 


INT29560 


C 


WRITE (6 , 300) ( YO( I) , 1=1, Nl) 


INT29570 


C 




INT29580 




RETURN 


INT29590 


C - • 




- - INT29600 


100 


F0RMAT( 715) 


INT29610 


200 


F0RMAT(6F10. 0) 


INT29620 


290 


F0RMAT( / , ' ORIGINAL COORDINATES ' , / , ' X/C ' ) 


INT29630 


295 


FORMAT ( / , ' INTERPLATED COORDINATES ' , / , ' X/C ' ) 


INT29640 


298 


FORMAT ( ’ Y/C ' ) 


INT29650 


300 


FORMAT( 6F10. 6) 


INT29660 


C - ■ 




- - INT29670 




END 


INT29680 


C 
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SUBROUTINE STAGR( N , STAG , XO , YO , XSTGR , YSTGR) 


INT29700 


C 




INT29710 




DIMENSION X0( 100) ,YO( 100) ,XSTGR( 100) , YSTGR( 100) ,DS( 100) 


INT29720 


C 




INT29730 




XOTE = 0.5 * (X0( 1)+X0(N) ) 


INT29740 




YOTE = 0.5 * (Y0(1)+Y0(N)) 


INT29750 




DS( 1) = SQRT((X0(1)-X0TE)**2 + ( Y0( 1) -Y0TE)**2) 


INT29760 




DSM = DS( 1) 


INT29770 




DO 10 I = 2 , N 


INT29780 




DS( I ) = SQRT( (X0( I ) -X0TE)**2 + ( Y0( I ) -YOTE) **2) 


INT29790 




IF (DS( I) .LT. DSM) GOTO 10 


INT29800 




IM = I 


INT29810 




DSM = DS( I ) 


INT29820 


10 


CONTINUE 


INT29830 


C 




INT29840 




YYY = YOTE-YO(IM) 


INT29850 




XXX = XOTE-XO(IM) 


INT29860 




IF (YYY .EQ. 0.0 .AND. XXX . EQ. 0.0) THEN 


INT29870 




ANG = 0. 0 


INT29880 




ELSE 


INT29890 




ANG = A T AN 2 ( Y” ', 7 , XXX ) 


INT29900 




END IF 


INT299 10 




ANG = ANG + STAG 


INT29920 


C 




INT29930 




COS AN = COS (ANG) 


INT29940 




SI NAN = SIN(ANG) 


INT29950 




DO 20 I = 1 , N 


INT29960 


c 


YY = Y0( I) -Y0( IM) 


INT299 7 0 


c 


XX = XO(I)-XO(IM) 


INT29980 


c 


IF (YY .EQ. 0.0 .AND. XX . EQ. 0.0) THEN 


INT29990 


c 


ANGCO = 0.0 


INT30000 


c 


ELSE 


INT30010 


c 


ANGCO = ATAN2( YY,XX) 


INT30020 


c 


END IF 


INT30030 




XSTGR( I )= X0( I )*COSAN + YO(I)*SINAN 


INT30040 




YSTGR( I )= Y0( I )*COSAN - XO(I)*SINAN 


INT30050 


20 


CONTINUE 


INT30060 



RETURN 

END 



INT30070 

INT30080 



APPENDIX B. C4 CASCADE 



A. EXPERIMENTAL RESULTS 

The experimental results of the C4 cascade were obtained directly from professor 
G.J. Walker, University of Tasmania, Tasmania, Australia, who performed these exper- 
iments. 

The results of the boundary layer measurements of the C4 cascade are given below 
at four inlet angles: 34.1°, 36.3°, 45.6°, and 47.7°. The Reynold numbers, based on the 
chord and the upstream velocity, are 200000, 191000, 173000 and 171000 respectively. 
The results given in the following tables include the displacement thickness (<5’), the 
shape factor (II) and the local free stream velocity (L'E). 



Table 1. EXPERIMENTAL RESULTS AT INLET ANGLE OF 34.1° 



x'c 


4* [10~ 3 FT] 


H 


UE [FT/SEC] 


0.4 


4.9 


2.48 


168.37 


0.5 


6.28 


2.61 


167.35 


0.6 


8.79 


3.24 


158.31 


0.7 


10.83 


3.63 


149.27 


O.S 


16.63 


3.79 


147.13 


0.9 


16.19 


1.89 


143.79 



Table 2. EXPERIMENTAL RESULTS AT INLET ANGLE OF 36.3° 



X c 


(5*|10“ J FT] 


H 


UE [FT/SEC] 


0.4 


5.43 


2.55 


161.63 


0.5 


7.09 


2.70 


157.59 


0.6 


10.3 


3.34 


148.78 


0.7 


12.63 


3.78 


139.87 


O.S 


14.84 


2.78 


135.01 


0.9 


16.43 


1.76 


133.23 
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Table 3. EXPERIMENTAL RESULTS AT INLET ANGLE OF 45.6° 



X c 


d [ 1 0 -3 FT ] 


H 


UE (FT, SEC] 


0.4 


8.08 


2.58 


137.88 


0.5 


9.83 


2.41 • 


133.70 


0.6 


12.35 


2.33 


122.18 


0.7 


12.98 


1.97 


114.93 


0.8 


19.44 


1.90 


111.77 


0.9 


27.69 


1.92 


109.26 



Table 4. EXPERIMENTAL RESULTS AT INLET ANGLE OF 47.7° 



X c 


5* [10~ 4 FT] 


H 


UE [FT/SEC] 


0.4 


8.87 


2.24 


130.64 


0.5 


10.27 


2.19 


124.76 


0.6 


14.31 


2.08 


116.86 


0.7 


16.45 


1.S7 


106.72 


0.8 


24.16 


1.82 


103.75 


0.9 


36.00 


2.01 


102.18 
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The results of the measurements of the velocity profiles in the boundary layer at two 
inlet angles. 34.1° and 36.3° at 50% chord are given below. 



Table 5. VELOCITY PROFILES AT 50% CHORD. 



y 


p = 36.3° 


p = 34.1° 


0.0 


0.0 


0.0 


2.3 


0.172 


0.208 


3.7 


0.270 


0.327 


6.2 


0.469 


0.534 


8.6 


0.666 


0.728 


11.0 


0.794 


0.867 


13.4 


0.891 


0.933 


18.3 


0.982 


0.985 


23,2 


1.000 


1 .000 



B. C4 CASCADE COORDINATES 



DIMENSION X( 0 : 100),XU(0: 100),XL(0: 100),YU(0: 100),YL(0: 100) C4 

DATA Al, A2, A3, A4/0. 15492,0. 06563,0. 2528,0. 2811/ C4 

DATA Bl,B2,B3,B4/0. 03866,0. 07871,0. 1467,0. 03448/ C4 

PI = AC0SC-1. 0) C4 

C READ (5,800) NMAX C4 

800 FORMAT (15) C4 

NHAX=33 C4 

C READ (5,810) (X( I) , 1=0, NMAX) C4 

810 FORMAT (6F10.6) C4 

DO 50 1=0, NMAX C4 

X(I) = (1. 0-COS(PI*I/NMAX))/2. C4 

50 CONTINUE C4 

DO 100 1=0, NMAX C4 

SRT = SQRT( (0. 5/SIN(PI/12))**2-(0. 5-X( I ) )**2) C4 

YC = -0. 5/TAN(PI/12) + SRT C4 

DY = ATAN( ( 0. 5 -X( I) )/SRT) C4 

IF (X( I). LT. 0. 3) THEN C4 

YT = A1' V SQRT(X( I) ) - A2*X(I) - A3*X(I)**2 + A4*X(I)**3 C4 

ELSE C4 

YT = B1 + B2*X( I) - B3*X(I)**2 + B4*X(I)**3 C4 

END IF C4 



00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 



o o 



YU(I) = YC + COS(DY)*YT C4 00220 

YL( I ) = YC - COS(DY)*YT C4 00230 

XU(I) = X(I) - SIN(DY)*YT C4 00240 

XL( I ) = X(I) + SIN(DY)*YT C4 00250 

100 CONTINUE C4 00260 

WRITE (6,900) ( I ,X( I) ,XU( I ) , YU( I ) ,XL( I) , YL( I ) , 1=0 ,NMAX) C4 00270 

900 FORMAT ( 15 , 4X, F10. 6 , 4X, 2F10. 6 , 4X, 2F10. 6) C4 00280 

WRITE (1,910) (XL( I) , I=NMAX,0 , -1),(XU(I), 1=1 ,NMAX) C4 00290 

WRITE (1,910) ( YL( I ) , I=NMAX, 0 , -1) , ( YU( I ) , 1=1 ,NMAX) C4 00300 

910 FORMAT (6F10.6) C4 00310 

STOP C4 00320 

END C4 00330 
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